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Abstract 


A computer  program  has  been  written  to  facilitate  performance  assessment  of  concrete  vaults 
used  in  Low  Level  Waste  (LLW)  disposal  facilities.  The  computer  program  is  a numerical 
computer  model  of  degradation  in  concrete.  A one-dimensional  finite  difference  equation 
is  used  to  propagate  ions  by  precipitation/dissolution  of  available  salts.  The  precipita- 
tion/dissolution of  salts,  in  turn,  changes  the  transport  properties,  which  changes  the  rate 
of  ion  transport.  The  result  is  a model  which  incorporates  the  synergism  of  multiple  degra- 
dation mechanisms. 

This  Report  is  self-contained.  It  includes  the  installation  instructions,  user  manual, 
technical  details,  and  source  code.  The  program  was  written  using  a literate  programming 
tool  and  the  “pretty-printing”  output  of  the  source  code  is  attached  at  the  end  of  this  report. 

Keywords:  building  technology;  computer  modelling;  concrete;  corrosion  of  reinforcement; 

degradation;  leaching;  radioactive  waste;  service  life;  sulfate  attack 
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1 Introduction 


The  computer  program  4SIGHT  has  been  written  to  facilitate  Low  Level  Waste  (LLW)  dis- 
posal facility  performance  assessment.  4SIGHT  can  be  used  to  predict  the  service  life  and 
hydraulic  conductivity  of  a buried  concrete  vault.  This  estimate  is  based  upon  an  analysis 
of  the  concrete  vault  roof  given  the  material  properties  of  the  concrete  and  the  external  ion 
concentrations. 

4SIGHT  incorporates  multiple  degradation  mechanisms  by  using  a single  transport  equa- 
tion for  ions.  Individual  degradation  mechanisms  are  typically  controlled  by  the  concentra- 
tion of  a single  ion  species.  However,  as  various  ion  species  diffuse  into  the  concrete,  the 
transport  properties  change,  which  changes  the  rate  of  ion  diffusion.  As  the  pore  solution 
pH  changes,  available  salts  either  precipitate  or  go  into  solution.  The  precipitation  and 
dissolution  processes  changes  the  porosity,  which  in  turn  changes  the  transport  properties. 
It  is  the  act  of  maintaining  chemical  equilibrium  that  affects  the  synergism  of  the  different 
degradation  mechanisms. 

Previous  efforts  [1] [2]  have  yielded  mathematical  models  for  various  degradation  mech- 
anisms. From  these  equations,  the  service  life  could  be  calculated  for  any  one  degradation 
mechanism.  However,  combining  the  effects  of  multiple  mechanisms  was  not  feasible  ana- 
lytically. The  advantage  of  the  method  used  by  4SIGHT  is  the  ability  to  incorporate  the 
synergism  of  multiple  degradation  processes.  Rather  than  calculate  rates  of  attack,  4SIGHT 
uses  a single  transport  equation  to  propagate  ions  through  the  concrete.  Keeping  the  system 
in  chemical  equilibrium  allows  interaction  of  degradation  mechanisms. 

4SIGHT  is  a command-driven  program  that  awaits  input  parameters  from  the  user,  using 
default  values  when  necessary.  Once  all  known  parameters  are  entered,  4SIGHT  reiterates 
the  values  of  all  the  input  parameters  used  in  the  calculation  with  a note  as  to  whether  the 
value  was  specified  by  the  user  or  whether  a default  value  was  used.  The  computations  begin 
by  4SIGHT  propagating  ions  in  discrete  time  increments.  At  regular  intervals  4SIGHT  prints 
to  the  screen  the  bulk  hydraulic  conductivity  and  diffusivity,  and  the  depth  of  sulfate  and 
chloride  attack.  The  calculation  continues  until  one  of  the  following  occurs:  corrosion  of  the 
steel  reinforcement,  sufficient  sulfate  attack  that  the  roof  fails  structurally,  failure  of  a joint, 
or  a user  specified  time  limit  has  been  exceeded.  Upon  termination,  4SIGHT  prints  the  state 
of  the  system,  composed  of  the  concentration  of  user  selected  ions  and  moles  of  user  selected 
salts,  as  a function  of  depth  into  the  roof. 

What  follows  in  this  Report  is  a complete  user’s  guide  to  4SIGHT,  including  installation 
information  and  source  code.  The  goal  is  an  all-inclusive  document,  acting  as  a single  body 
of  work,  that  facilitates  review  of  any  aspect  of  the  program.  To  facilitate  source  code 
review,  4SIGHT  was  written  using  a literate  programming  tool  which  incorporates  both  the 
programming  code  and  typeset  comments.  The  typeset  output  does  not  conform  to  the 
structure  of  a NISTIR.  Therefore,  the  typeset  output,  as  it  normally  appears,  is  attached  to 
the  end  of  this  report. 
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2 Installation 


2.1  Distributed  Files 

4SIGHT  is  distributed  with  the  following  files: 


README 

4SIGHT.EXE 

ION.DB 

MANUAL.PS 

MANUAL . HP 

MANUAL.TXT 

EXAMPLE.DAT 


Basic  installation  instructions 
Executable  program 
Ion  database 

4SIGHT  manual  in  PostScript  format 
4SIGHT  manual  in  HP-PCL  format 
4SIGHT  manual  in  ASCII  format 
Example  input  file 


To  install  4SIGHT,  simply  copy  these  files  to  a directory  on  the  hard  drive.  This  manual  has 
been  included  in  three  different  useful  formats  for  printing. 


2.2  Hardware  Requirements 

To  execute  the  program  4SIGHT,  the  user  must  be  using  a personal  computer  equipped  with 
either  an  80386  or  80486  microprocessor  with  a math  coprocessor.  Lacking  either  of  these 
requirements,  the  program  will  fail  to  run. 


3 Program  Execution 

4SIGHT  was  written  to  be  executed  in  either  interactive  or  batch  mode.  The  interactive 
mode  simply  waits  for  the  user  to  type  in  the  various  input  parameters  at  run-time.  When 
executed  in  batch  mode,  the  user  first  edits  an  ASCII  file  containing  the  input  parameters 
as  if  typed  during  interactive  mode.  This  ASCII  file  is  then  included  at  the  command  line 
when  running  4SIGHT. 

3.1  Interactive  Mode 

When  run  interactively,  4SIGHT  simply  waits  for  the  user  to  enter  input  parameters.  The 
user  specifies  the  end  of  the  list  by  entering  either  quit  or  exit.  For  example,  assuming  the 
executable  program  is  in  the  directory  C:\4SIGHT,  interactive  mode  is  initiated  by  simply 
typing  the  program  name: 


C:\4SIGHT\4SIGHT 


The  program  responds  with: 
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Enter  commands: 


At  this  point  4SIGHT  is  simply  waiting  for  the  user  to  input  program  parameters.  The 
commands  are  given  one  per  line: 

DIFF  = 2.0E-12 
WC  = 0.45 
TIME  = 100000 
quit 

After  which,  4SIGHT  begins  the  calculation,  printing  intermediate  results  to  the  screen. 

3.2  Batch  Mode 

When  run  in  batch  mode  the  user  specifies  an  input  file  at  the  command  line.  This  input  file 
simply  contains  any  number  of  input  parameters,  as  in  the  interactive  mode.  For  example, 
assume  the  following  text  is  saved  into  a file  called  input  .dat: 

DIFF  = 5.0E-12 
THICKNESS  =1.0 
WC  = 0.50 

EXTERNAL  Cl  = 0.150 
EXTERNAL  Na  = 0.150 
DEPTH  =0.25 

This  file  can  be  created  by  any  ASCII  editor  such  as  the  DOS  edit  program  or  the  UNIX 
vi  program.  When  using  a commercial  word  processing  program  to  create  input  files  make 
sure  that  the  output  is  in  plain  ASCII.  As  a check,  give  the  C:\type  input.dat  command 
at  the  DOS  prompt.  The  text  should  appear  as  typed.  Note  that  neither  quit  nor  exit  are 
required  in  the  input  file.  To  use  input  files  with  4SIGHT,  simply  include  that  file  name  on 
the  DOS  command  line: 

C:\4SIGHT\4SIGHT  input.dat 

The  output  from  4SIGHT  goes  to  the  screen.  This  is  useful  for  initial  calculations  and  trial 
runs.  However,  if  the  user  desires  to  save  the  output  to  an  ASCII  file,  simply  use  redirection 
of  the  output: 

C:\4SIGHT\4SIGHT  input.dat  >output.dat 

This  example  uses  the  commands  in  the  file  input . dat  as  input  parameters  and  saves  the 
output  in  the  file  output  .dat. 
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3.3  Plotting  Results 

The  output  from  4SIGHT  has  been  formatted  with  quotes  (")  and  tabs  to  facilitate  incorpo- 
rating them  into  commercial  spreadsheet  programs.  Once  incorporated  into  the  spreadsheet 
program,  the  user  can  plot  either  the  sulfate  and  chloride  depth  versus  time  from  the  inter- 
mediate results,  or  the  ion  concentration  as  a function  of  depth  from  the  final  state  of  the 
system. 


4 Input  Specifications 

The  complete  list  of  input  parameters  is  given  here  to  introduce  the  parameters  used  in  the 
example  to  follow.  Although  many  of  the  parameters  are  self-evident,  the  complete  meaning 
of  all  the  parameters  will  not  be  understood  until  the  user  reads  the  explanatory  sections 
that  follow  the  example. 

The  input  commands  to  4SIGHT  are  straightforward.  With  very  few  exceptions,  the 
format  of  the  commands  is 


parameter  = value 

Even  though  every  possible  input  parameter  is  listed  below,  there  are  sufficient  default 
parameters  for  the  user  to  perform  a calculation  without  entering  any  parameters. 

The  syntax  of  the  list  of  parameters  below  is  as  follows:  Text  to  be  typed  verbatim  is 
shown  in  typewriter  font  ( e.g .,  DIFF  = ).  These  words  can  be  typed  in  any  mixture  of  upper 
and  lower  case  characters  with  any  number  of  spaces  or  tabs  between  the  parameters  and 
the  ‘=\  Words  in  italics  represent  user  selected  input  values  (e.g.,  expr).  There  are  two  data 
types:  numbers  and  ions,  which  are  represented  by  expr  and  ion , respectively.  Numbers 
can  be  given  in  integer,  fixed  point,  or  scientific  notation.  Examples  of  valid  values  for  expr 
include: 


120  120.0  1.2E+02 

The  units  for  each  input  value  appear  under  expr  in  the  definition.  4SIGHT  also  recognizes 
the  following  ions  (the  valence  has  been  omitted): 

H Ca  Na  K OH  Cl  S04  C03 

The  ion  values  MUST  appear  as  shown  since  ion  variables  are  case  sensitive.  4SIGHT  can 
also  recognize  salts.  A salt  is  specified  by  concatenating  a cation  and  an  anion,  separated 
by  a space.  Stoichiometric  ratios  are  not  used.  For  example,  the  following  specifications 

Na  Cl  Ca  OH  Na  S04 
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are  for  sodium  chloride,  NaCl , calcium  hydroxide,  Ca(OH) 2,  and  sodium  sulfate,  Na2S04 , 
respectively. 

Below  is  a list  of  all  possible  input  parameters.  At  the  end  of  each  description  is  an 
example  of  usage.  In  cases  where  a default  value  exists,  the  default  value  is  used  as  the 
example  value.  Whether  or  not  the  value  is  a default  value  is  noted  next  to  the  example. 


4.1  Material  Properties 


DIFF  = expr 

m2 1 sec 


expr : 


default: 


The  ultimate  Cl~  diffusivity  of  the  concrete.  This  value  represents  the 
diffusivity  of  the  undamaged  concrete.  The  default  value  given  below  is 
suitable  for  a w/c=0.45  ordinary  portland  concrete. 

DIFF  = 5.7E-12 


PERM  = expr 

,2 


m 


expr: 


default: 


The  ultimate  permeability  of  the  concrete, 
permeability  of  the  undamaged  concrete. 
PERM  = 2.5E-18 


This  value  represents  the 


WC 


expr 


expr:  Concrete  water:cement  weight  ratio  during  mixing. 

default:  WC  = 0.45 

4.2  Ion/Fluid  Initial  Conditions 

EXTERNAL  ion  = expr 
mol/L 

ion:  The  ion  of  interest. 

expr:  The  external  concentration  of  ion  above  the  roof.  This  command  can 

be  repeated  for  each  ion  external  to  the  roof  slab. 
example:  EXTERNAL  Cl  = 0.15 

INTERNAL  ion  = expr 
mol/L 
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ion : 

expr: 

Ion  of  interest. 

The  concentration  of  ion  in  the  pore  solution  of  the  concrete.  This 
command  was  designed  to  let  the  user  specify  the  internal  Na+  and 
K+  concentration  to  establish  the  correct  pH.  This  command  can  be 

example: 

repeated  for  each  relevant  ion. 

INTERNAL  Na  = 0.10 

HEAD  = expr 
m 


expr: 

The  equivalent  height  of  ground  water  (hydraulic  head)  above  the  upper 

default: 

surface  of  the  roof. 

HEAD  =5.00 

4.3  Geometry. 

THICKNESS  = expr 
m 

expr : The  vertical  thickness  of  the  roof. 

default : THICKNESS  =1.00 

CRACK  = exprl  AT  expr2  DEPTH  exprS 
mm  m 


exprl : 
expr 2: 
expr 3: 
default: 

The  width  of  the  crack. 

The  spacing  between  cracks. 

The  penetration  of  the  crack,  relative  to  the  bottom  of  the  slab. 

CRACK  = 0.000100  AT  1.00  DEPTH  0.30 

JOINT  PERM  = 

expr 

9 

m 

expr: 

The  permeability  of  the  joint  filling  compound  between  roof  slabs.  This 
value  should  not  reflect  the  geometry  of  joint.  Rather,  it  is  a material 

example: 

property  of  the  joint  compound. 

JOINT  PERM  = 3.2E-15 

JOINT  = exprl  AT  expr2  UNTIL  exprS 
m m years 
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exprl : 

The  width  of  the  joint. 

expr 2: 
expr 3: 
example: 

The  spacing  between  joints. 

The  service  life  of  the  joint  compound. 

JOINT  = 0.020  AT  1.00  UNTIL  100 

4.4  Failure/Termination  Limits 

REBAR  = expr 
m 


expr: 

The  depth  of  the  reinforcement  bars  (rebars)  from  the  top  of  the  slab. 

example: 

When  a critical  concentration  of  Cl~  ions  has  penetrated  to  a depth 
expr  the  structure  fails. 

REBAR  =0.90 

DEPTH  = expr 


m 

expr: 

The  critical  penetration  depth  for  sulfate  attack.  Once  the  sulfate  front 
has  penetrated  down  to  a depth  expr  from  the  top,  the  roof  is  unable  to 

default: 

sustain  its  load  and  it  fails. 

DEPTH  =0.20 

TIME  = expr 
day 

expr: 

The  maximum  number  of  days  to  continue  the  calculation.  Note  that 

default: 

no  comma  should  be  used. 

TIME  = 100000 

4.5  Sulfate  Attack  Parameters 

YOUNGS  = expr 
GPa 

expr:  The  Youngs  modulus  of  the  concrete. 


default: 

YOUNGS  =20.0 

BETA  = expr 


expr: 

default: 

The  linear  strain  due  to  one  mole  of  Na+  per  m3  of  concrete. 

BETA  = 1.8E-06 
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CE  = expr 

mol/m3 

expr : The  concentration  of  reacted  sulfate  as  ettringite. 


default: 

CE  = 350.0 

ROUGHNESS  = 

expr 

expr: 

default: 

The  fracture  surface  roughness  factor. 

ROUGHNESS  =1.0 

GAMMA  = expr 
J/m2 

expr:  The  fracture  surface  energy  of  the  concrete. 

default:  GAMMA  =10.0 

POISSON  = expr 


expr: 

The  Poisson  ratio  of  the  concrete. 

default: 

POISSON  = 0.2 

4.6  Output  Parameters 

OUTPUT  = ion 


ion: 

Output  the  pore  solution  concentration  ( mol/L ) of  ion  in  the  output 
describing  the  final  state  of  the  system.  This  command  lets  the  user 
examine  the  concentration  of  specified  ions  as  a function  of  depth.  This 

example: 

command  can  be  repeated  for  different  ions. 

OUTPUT  = C03 

OUTPUT  = ionl  ion2 


NOTE: 

Output  the  quantity  (moles)  of  salt  in  the  volume  Vsample.  The  salt  is 
comprised  of  the  cation  ionl  and  the  anion  ion2. 

The  two  ion  specifications  are  separated  by  a space  and  stoichiometric 

example: 

ratios  are  omitted. 

OUTPUT  = Ca  OH 
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5 Example 

To  illustrate  the  use  of  4SIGHT,  consider  the  following  fictitious  example.  The  concrete  was 
designed  with  a 0.40  waterxement  ratio  (WC=0.40).  Experimental  diffusivity  measurements 
using  Cl~  ions  gave  5.0  x 10 ~12  m2 / sec  (DIFF=5 . OE-12).  No  permeability  measurements 
are  available.  The  roof  slab  is  1.0  meter  thick  (THICKNESS=1 . 000).  Regularly  spaced  cracks 
have  been  observed  on  the  bottom  of  the  slab.  The  cracks  are  approximately  100  [im  wide, 
spaced  2 meters  apart,  and  are  assumed  to  extend  upwards  to  the  neutral  axis  which  is 
25  centimeters  from  the  bottom  of  the  roof  slab  (CRACK  = 0.000100  AT  2.0  DEPTH  0.25). 
Also,  the  roof  is  buried,  giving  an  effective  pressure  head  of  2 meters  (HEAD  = 2.0). 

Soil  analysis  indicates  the  presence  of  SO\~  at  a concentration  of  1.0  moles  per  liter 
(EXTERNAL  S04  = 1.0).  Engineering  analysis  indicates  that  if  the  sulfate  degradation  pen- 
etrates down  20  centimeters  from  the  top  surface  of  the  roof  then  the  vault  will  collapse 
(DEPTH  = 0.20).  Additionally,  chloride  ions  are  present  in  the  soil  at  a concentration  of  0.40 
moles  per  liter  (EXTERNAL  Cl  = 0 . 40)  and  engineering  drawings  indicate  that  the  reinforce- 
ment bars  are  located  60  centimeters  from  the  top  of  the  roof  (REBAR  = 0.60).  Thorough 
soil  analysis  indicates  that  sodium  ions  are  also  present  at  a concentration  of  2.40  moles  per 
liter  (EXTERNAL  Na  = 2.40),  giving  a nearly  neutral  soil  pH. 

Internal  to  the  concrete  the  pH  is  approximately  13.  Therefore,  internal  potassium  and 
sodium  concentrations  are  approximately  0.1  (INTERNAL  K = 0.1)  and  0.05  (INTERNAL  Na 
= 0.05)  moles  per  liter,  respectively. 

To  monitor  the  ingress  of  the  external  ions,  calcium  (OUTPUT  Ca)  and  chloride  (OUTPUT 
Cl)  ions  will  be  included  in  the  output  of  the  final  state  of  the  system.  Also,  to  monitor 
leaching,  the  solid  calcium  hydroxide  content  will  also  be  included  in  the  output  (OUTPUT  Ca 
OH). 

The  input  file  for  this  example,  example.dat  is  included  in  the  distribution  diskette  and 
reiterated  here: 

DIFF  = 5. OE-12 
WC  = .40 

THICKNESS  = 1.000 
EXTERNAL  Na  = 2.40 
EXTERNAL  Cl  = 0.400 
EXTERNAL  S04  =1.00 
INTERNAL  K = 0.1000 
INTERNAL  Na=  0.0500 
OUTPUT  Ca 
OUTPUT  Cl 
OUTPUT  Ca  OH 
REBAR  = .8000 
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HEAD  =2.0 

CRACK  = 0.000100  AT  2.0  DEPTH  0.25 
DEPTH  = .2000 
TIME  = 100000 


To  use  this  input  file,  simply  include  example.dat  at  the  DOS  command  line.  To  save  the 
results,  redirect  the  output  to  a file: 

C:\4SIGHT\4SIGHT  EXAMPLE.DAT  >EXAMPLE.0UT 

In  this  example,  the  output  is  stored  in  EXAMPLE. OUT,  an  ASCII  file  which  can  be  imported 
into  any  spreadsheet  program. 

NOTE:  This  calculation  required  90  seconds  to  complete  on  a personal  computer  equipped 
with  a 80486  microprocessor  operating  at  66MHz. 


5.1  Internal  Parameters. 

4SIGHT  outputs  all  of  the  parameters  the  user  has  specified  ("USER"),  along  with  default  and 
calculated  values  ("DEFAULT"). 

" This  is  4SIGHT  (Version  1.0)" 


"THICKNESS 

II 

1.00000 

"(m) 

" "USER" 

"DIFF 

II 

5 . Oe-12 

" (m~2/ sec) 

" "USER" 

"PERM 

II 

9 . 8e-13 

" (m/ sec) 

" "DEFAULT" 

"WC 

II 

0.40000 

II 

" "USER" 

"HEAD 

II 

2.00000 

"(m) 

" "USER" 

Sulfate  Attack 

Parameters : 

"YOUNGS 

II 

2 . 0e+10 

" (N/nT2) 

" "DEFAULT" 

"BETA 

II 

1 . 8e-06 

II 

" "DEFAULT" 

"CE 

II 

350.0000C 

"(Mol/m‘3)"  "DEFAULT 

"ROUGHNESS 

II 

1.00000 

II 

" "DEFAULT" 

"GAMMA 

II 

10.00000 

" ( J/m~2) 

" "DEFAULT" 

"POISSON 

II 

0.20000 

II 

" "DEFAULT" 

"DEPTH 

II 

0.20000 

" (m) 

" "USER" 

"REBAR 

II 

0.80000 

" (m) 

" "USER" 

"TIME 

II 

100000 

" (day) 

" "USER" 

"CRACK  = 0.00010  AT  2.00000  DEPTH  0.25000" 


"Chloride  failure  (yr)"  1743 

"Sulfate  failure  (yr)"  447 
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. 


5.2  Initial  State 


The  report  of  the  initial  state  lets  the  user  verify  the  EXTERNAL  and  INTERNAL  conditions  of 
the  system. 


Initial 

state  of 

system : 

II 

ION 

II 

"EXTERNAL 

" "INTER 

II 

H 

II 

0.00000 

0.00000 

II 

Ca 

II 

0.00000 

0.00035 

II 

Na 

II 

2.40000 

0.05000 

II 

K 

II 

0.00000 

0.10000 

II 

OH 

II 

0.00000 

0.15070 

II 

Cl 

II 

0.40000 

0.00000 

II 

S04 

II 

1.00000 

0.00000 

II 

pH 

II 

7.00000 

13.17810 

These  results  indicate  that  the  pH  of  the  environment  is  7,  while  the  pH  of  the  pore 
solution  is  initially  13.  To  change  the  pH  of  the  environment  the  user  can  simply  change  the 
concentration  of  EXTERNAL  anions  or  cations. 

5.3  Depth  vs.  Time 

As  4SIGHT  is  calculating  ion  transport  it  regularly  prints  the  current  status  of  the  degra- 
dation. The  sulfate  and  chloride  penetration  depths  as  a function  of  time  are  given  in  the 
columns  labelled  S04  and  Cl,  respectively.  L is  the  remaining  thickness  of  the  slab,  K is  the 
hydraulic  conductivity,  D is  the  diffusivity,  Flux  is  the  flux  of  pore  solution  out  the  bottom 
of  the  slab,  and  pH  is  the  pH  of  the  flux. 
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"Day" 

"L" 

"K  " 

"D 

" "S04" 

"Cl" 

"Flux"  "pH" 

II  II 

"m" 

"m/s" 

"m~2/s 

" "m  " 

"m  " 

"ml/dy/m2" 

0 

1.00C 

i 1 . 3e-12  5 . Oe-12 

0.000 

0.000 

0.027 

13.2 

10036 

0.988 

1 .3e-12 

5. Oe-12 

0.012 

0.101 

0.027 

13.2  - 

15055 

0.982 

1.3e-12 

5. Oe-12 

0.018 

0.129 

0.027 

13.2 

20023 

0.976 

1 .3e-12 

5. Oe-12 

0.025 

0.150 

0.028 

13.2 

25031 

0.969 

1 . 3e-12 

5. Oe-12 

0.031 

0.173 

0.028 

13.2 

30018 

0.963 

1 . 3e-12 

5. Oe-12 

0.037 

0.192 

0.028 

13.2 

35011 

0.957 

1.3e-12 

5. Oe-12 

0.043 

0.211 

0.028 

13.2 

40048 

0.951 

1 . 3e-12 

5. Oe-12 

0.049 

0.229 

0.029 

13.2 

45104 

0.945 

1 . 3e-12 

5. Oe-12 

0.055 

0.246 

0.029 

13.2 

50017 

0.939 

1 ,3e-12 

5. Oe-12 

0.061 

0.263 

0.029 

13.2 

55045 

0.933 

1 . 3e-12 

5. Oe-12 

0.067 

0.280 

0.029 

13.2 

60015 

0.927 

1 . 3e-12 

5. Oe-12 

0.074 

0.295 

0.030 

13.2 

65003 

0.920 

1 . 3e-12 

5. Oe-12 

0.080 

0.311 

0.030 

13.2 

70002 

0.914 

1 . 3e-12 

5. Oe-12 

0.086 

0.327 

0.030 

13.2 

75013 

0.908 

1 . 3e-12 

5. Oe-12 

0.092 

0.342 

0.031 

13.2 

80069 

0.902 

1.4e-12 

5. Oe-12 

0.098 

0.357 

0.031 

13.2 

85037 

0.896 

1.4e-12 

5. Oe-12 

0.104 

0.372 

0.031 

13.2 

90050 

0.890 

1 . 4e-12 

5. Oe-12 

0.110 

0.387 

0.031 

13.2 

95069 

0.884 

1 . 4e-12 

5. Oe-12 

0.116 

0.401 

0.031 

13.2 

100017 

0.878 

l 1.4e-12  5. Oe-12 

0.122 

0.416 

0.032 

13.2 

As  per  the  linear  model  for  sulfate  attack,  the  sulfate  front  increases  linearly  with  time. 
The  chloride  depth  has  a t1^2  because  of  the  low  Peclet  number  (the  ratio  of  Darcy  to 
‘diffusive’  flow  [2]).  At  sufficiently  high  Peclet  numbers  (greater  permeability  or  hydraulic 
head)  the  chloride  depth  approaches  a linear  relationship  to  time. 


5.4  Failure  Data 

After  reporting  the  time  dependent  behavior,  4SIGHT  reports  the  reason  for  termination  and 
the  status  of  the  sulfate  and  chloride  penetration: 

Exceeded  TIME  limit. 

"T"  17.283 

"Day"  100017 
"S04  (m)"  0.122 

"Cl  (m)"  0.416 
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These  failure  data  indicate  that  the  calculation  terminated  because  it  exceeded  the  time 
limit.  The  termination  occurred  after  100044  days  (586  years),  at  which  point  the  sulfate 
and  chloride  penetration  were  0.13  and  0.60  meters,  respectively. 


5.5  Final  System  State 

Upon  termination  of  the  calculation,  the  final  state  of  the  system  for  the  slab  is  printed.  The 
L(m)  is  measured  from  the  top  of  the  slab.  Psi  is  the  hydraulic  pressure,  vD  is  the  Darcy 
velocity,  xi  is  the  inverse  of  the  formation  factor,  phi  is  the  porosity,  f c is  the  estimated 
compressive  strength  using  ACI  211. 


Final  System 
"Psi" 


"L(m) " 
0.0000 
0.0500 
0.1000 
0.1500 
0.2000 
0.2500 
0.3000 
0.3500 
0.4000 
0.4500 
0.5000 
0.5500 
0.6000 
0.6500 
0.7000 
0.7500 
0.8000 
0.8500 
0.9000 
0.9500 
1.0000 


0.392 

0.392 

0.392 

0.374 

0.342 

0.311 

0.279 

0.248 

0.217 

0.186 

0.155 

0.124 

0.093 

0.062 

0.031 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 


state 

"vD" 

0.000 

0.000 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 

0.031 


"xi" 

0.0024 

0.0024 

0.0024 

0.0024 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 

0.0025 


"phi" 

0.1634 

0.1631 

0.1658 

0.1663 

0.1671 

0.1672 

0.1672 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 

0.1673 


"pH"  "fc"  "Ca" 
7.000  5052  0.0000 
7.000  5057  0.0000 
7.000  5011  0.0000 
13.277  5003  0.0001 
13.351  4989  0.0001 
13.392  4988  0.0001 
13.405  4987  0.0001 
13.394  4986  0.0001 
13.365  4986  0.0001 
13.326  4986  0.0002 
13.286  4986  0.0002 
13.250  4986  0.0002 
13.220  4986  0.0003 
13.199  4986  0.0003 
13.185  4986  0.0003 
13.176  4986  0.0004 
13.173  4986  0.0004 
13.172  4986  0.0004 
13.172  4986  0.0004 
13.172  4986  0.0004 
13.172  4986  0.0004 


"Cl"  "CaOH" 
0.4000  0.000 
0.4000  0.000 
0.4000  0.000 
0.3740  31.599 
0.3253  31.860 
0.2773  31.886 
0.2314  31.899 
0.1889  31.912 
0.1509  31.918 
0.1178  31.919 
0.0899  31.919 
0.0671  31.919 
0.0488  31.919 
0.0348  31.919 
0.0242  31.919 
0.0164  31.919 
0.0110  31.919 
0.0073  31.919 
0.0050  31.919 
0.0039  31.919 
0.0039  31.919 


Note  that  the  external  quantities  are  duplicated  to  a depth  of  0.1  meters.  As  the  sulfate 
penetrates  the  concrete,  the  concrete  fails  and  the  external  boundary  conditions  move  into 
the  concrete.  Also  note  that  the  pressure  (P)  is  zero  from  0.75  to  1.00  meters.  Since  the 
concrete  was  cracked,  the  permeability  of  the  cracked  portion  of  the  concrete  overwhelmed 
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the  uncracked  portion,  resulting  in  virtually  no  pressure  drop  across  the  crack.  This  increases 
the  pressure  gradient  across  the  remaining  uncracked  concrete. 


6 System 

Although  each  LLW  facility  may  be  unique,  most  underground  facilities  can  be  represented 
schematically  as  in  Figure  1.  Below  the  ground  surface  there  will  be  an  engineered  barrier 
to  deflect  surface  water  away  from  the  concrete  vault  below.  The  entire  facility  is  sited  on  a 
geologically  suitable  location. 


Figure  1:  Schematic  of  underground  LLW  facility 

From  the  schematic  in  Figure  1 it  appears  as  though  the  roof  is  the  most  critical  element 
because  it  is  most  likely  to  have  a moist  environment,  especially  upon  the  failure  of  the 
engineered  barrier.  Because  of  this,  the  analysis  of  the  entire  concrete  vault  can  concentrate 
on  the  roof.  If  the  roof  is  treated  as  a simple  slab,  the  analysis  simplifies  further.  Since  the 
flow  through  the  roof  slab  will  be  approximately  uniform  over  the  surface  of  the  slab,  a one- 
dimensional analysis  of  transport  vertically  through  the  slab  should  serve  as  a sufficiently 
accurate  model  for  transport  through  the  roof  of  the  vault. 

7 Ion  Transport 

A single  transport  equation  is  developed  to  propagate  ions  through  the  slab.  This  equation 
can  be  converted  into  a finite  difference  equation  so  that  it  can  be  implemented  in  a computer 
program.  A result  of  this  approach  is  that  time  will  advance  in  discrete  intervals.  After  every 
time  interval  each  computational  element  is  put  in  chemical  equilibrium  using  solubility 
products  and  a charge  balance.  This  step  is  achieved  through  dissolution/precipitation  of 
available  salts.  Any  change  in  the  quantity  of  solid  salts  in  the  pore  space  will  change  the 
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porosity  of  the  concrete,  hence  changing  the  transport  coefficients.  This  is  how  the  synergism 
of  degradation  mechanisms  is  achieved. 

7.1  Advection-Diffusion 

At  the  core  of  4SIGHT  is  the  advection-diffusion  equation  which  establishes  the  transport 
of  ions  through  the  slab.  The  advection-diffusion  equation  is  simply  the  diffusion  equation 
with  an  extra  term  to  account  for  Darcy  flow.  The  development  of  the  equation  starts  most 
simply  from  the  relation  between  flux,  j,  and  concentration,  c: 

j = —D'Ve  + cu  (1) 

The  parameter  D is  the  effective  diffusion  coefficient,  and  the  quantity  u is  the  average  pore 
fluid  velocity.  The  rate  of  change  in  concentration  is  the  negative  divergence  of  eqn.  1: 

dc 

— =V-DVc-u-Vc  (2) 

dt  v ’ 

after  neglecting  the  divergence  of  the  volume  averaged  velocity,  since  the  fluid  is  virtually 
incompressible  and  the  rate  change  in  porosity  is  small. 

7.2  Darcy  Flow 

The  average  pore  fluid  velocity  can  be  calculated  from  the  Darcy  velocity,  which  is  the 
volume- averaged  pore  fluid  velocity.  Given  a porous  media  with  permeability  k and  pore 
fluid  viscosity  p,  the  Darcy  velocity,  Vd,  is  proportional  to  the  pressure  gradient,  Vp  and 
the  density  of  the  fluid  [3]: 

Vd  = -- {Vp  - pg)  (3) 

The  Darcy  velocity  vp  can  be  related  to  the  average  pore  velocity,  u,  from  a geometrical 
argument.  Let  v(x)  represent  the  velocity  of  the  pore  fluid  at  some  point  x in  the  slab.  As- 
sume that  the  pore  fluid  completely  fills  the  available  pore  space.  Representing  the  porosity 
by  0,  u is  defined  as: 

u=ivlv{x)dV  (4) 

The  Darcy  velocity  is  the  average  over  the  entire  volume  V : 

vD  = i ^ v(x)  dR  (5) 

Since  v(x)  = 0 outside  of  the  porosity,  eqn.  5 can  be  limited  to  the  pore  space: 

vd  = -7  / v(x)  dV  (6) 

V J<t> 
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(7) 


From  comparison  to  eqn.  4,  the  relation  between  v#  and  u is: 

\D  = </>  u 


7.3  Continuity 

Once  the  Darcy  equation  is  incorporated  into  eqn.  2,  the  continuity  equation  is  needed 
in  order  to  update  the  pressures.  Upon  execution  of  the  computer  program,  the  follow- 
ing sequence  is  continuously  reiterated:  transport  ions,  attain  chemical  equilibrium,  adjust 
transport  properties,  update  boundary  conditions.  As  the  porosity  of  the  concrete  changes, 
the  transport  properties  also  change.  The  change  in  the  transport  properties  will  effect 
the  hydraulic  pressure  distribution  in  the  concrete  since  the  transport  properties  will  not 
by  uniform  throughout  the  concrete.  The  pressure  distribution  will  be  updated  using  the 
continuity  equation  for  porous  media. 

The  continuity  equation  for  a fluid  is 

% + V • pv  = 0 (8) 

To  develop  a continuity  equation  for  porous  media,  find  the  average  integral  value  of  eqn.  8 
over  a representative  volume,  V : 

vIvTtdV  + vIvv^vdV  = 0 (9> 

Rearranging  the  order  of  integration  gives 

UlvpdV  + ^-vIvPVdV  = 0 (10) 

Finally,  assuming  p is  constant  and  using  the  definition  of  Darcy  velocity  gives 

^ + V-vD  = 0 (11) 

This  is  the  result  obtained  by  Slattery  [4]  for  porous  media.  Substituting  for  vd  from  eqn.  3 
gives 

^7  = V • — (Vp  — pg)  (12) 

Ot  p 


7.4  Dimensionless  Variables 


Eqn.  2,  3,  and  11  can  be  combined  to  form  a system  of  equations  to  propagate  ions  through 
a porous  media.  However,  the  system  of  equations  can  be  condensed  by  a transformation 
into  dimensionless  variables.  Consider  the  following  definitions: 


P = 


V 

pD  o / k0 


(13) 
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where  the  capital  letters  A',  T,  and  P refer,  respectively,  to  the  dimensionless  distance, 
time,  and  pressure.  The  quantity  D0  is  the  chloride  diffusivity  of  the  concrete,  and  kQ  is  the 
permeability.  L is  simply  any  characteristic  length.  Using  these  definitions,  the  advection- 
diffusion  equation  (2)  becomes 


Svp+rc'Vc 


(14) 


This  equation  will  be  simplified  further  later  using  material  properties  related  to  porosity. 


7.5  Atkinson-Hearne  Model 


As  yet,  a dependable  mechanistic  model  does  not  exist  for  sulfate  attack.  Therefore,  to 
incorporate  sulfate  attack,  the  depth  of  sulfate  degradation  is  calculated  using  the  Atkinson 
and  Hearne  [5]  model: 


E(3C0CeD 
<27(1  — u) 


(15) 


E Youngs  modulus 

C0  External  sulfate  concentration 
CE  Concentration  of  sulfate  as  ettringite 
D Sulfate  diffusion  coefficient 
ot  Roughness  factor 

7 Fracture  surface  energy 

v Poisson  ratio 

Since  the  Atkinson-Hearne  model  gives  the  location  of  the  sulfate  front,  4SIGHT  presumes 
that  all  of  the  concrete  behind  the  sulfate  front  has  been  completely  disintegrated,  giving  it 
the  properties  of  the  surrounding  soil.  Since  the  transport  coefficients  of  soil  are  much  larger 
than  concrete,  the  external  boundary  conditions  are  advanced  to  the  sulfate  front,  creating 
a moving  boundary  condition. 

The  Youngs  modulus,  roughness  factor,  fracture  surface  energy,  and  Poisson  ratio  must  be 
determined  from  experimental  measurements.  The  quantity  CE  can  be  either  calculated  from 
the  cement  composition,  or  more  accurately  from  experiment.  To  experimentally  determine 
CE , the  sulfate  reacted  per  unit  mass  hydrated  cement  is  plotted  against  the  logarithm  of 
time.  This  data  is  fit  to  the  equation  [5] 


m = 


(16) 


where  m is  the  moles  of  sulfate  reacted  in  the  cement,  m0  is  the  free  parameter  of  the 
regression,  t is  time,  c is  the  concentration  of  sulfate  in  liquid,  tr  is  the  characteristic  time 
for  reaction,  and  is  the  concentration  in  kinetic  experiments.  The  maximum  value  of  m 
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can  be  calculated  from  the  initial  quantity  of  C3A  in  the  cement.  If  this  maximum  value 
is  labelled  mc,  the  value  of  t for  m = mc  is  the  time  at  which  all  of  the  C3A  is  consumed, 
t too • 

The  model  proposed  by  Atkinson  and  Hearne  assumes  that  reaction  is  the  controlling 
rate  process.  Therefore,  the  time  to  spalling  should  be  greater  than  the  time  required  for 
complete  consumption  of  the  C3A.  The  time  to  spalling  is  [5] 


where 


tspall  — 


Y2  — 

on/i  / / 


K^iCe 

2 Dc0 

(17) 

207(1  — u) 

E(PCe )2 

(18) 

In  the  event  that  tspaii  < ^ocn  Ce  must  be  calculated  in  a self  consistent  manner  using  eqns. 
16,17,  and  18.  More  complete  details  can  be  found  in  [5].  For  ordinary  portland  cements,  it 
is  likely  that  t3pau  > t^.  Therefore,  since  it  is  assumed  that  the  quantity  of  external  sulfate 
is  sufficiently  great  to  act  as  an  infinite  reservoir  (external  concentration  is  constant),  Ce 
can  be  calculated  from  the  C3A  content  of  the  cement.  Since  each  mole  of  ettringite  requires 
one  mole  of  A/203[6],  the  molar  concentration  of  AI2O3  will  be  the  molar  concentration  of 
ettringite.  Given  a concrete  mix  design  having  xcem  kilograms  of  cement  per  cubic  meter  of 
concrete,  and  the  cement  having  a weight  fraction  <f>Aho3  of  aluminum  oxide,  the  moles  of 
ettringite  formed  per  cubic  meter  of  concrete  is 


Ce 


%cem &AI2O3 

0.10196 


with  the  AI2O3  molar  weight  of  0.10196  kg  per  mole. 


(19) 


7.6  Effects  of  Cracks  and  Joints 

Although  not  a degradation  mechanism  per  se,  effects  due  to  the  presence  of  cracks  and 
joints  can  be  incorporated  into  the  models  for  transport  of  the  ions.  Since  the  roof  slab  will 
likely  be  a supported  member  such  that  the  bottom  of  the  roof  slab  will  be  in  tension,  it 
is  assumed  that  cracks  will  appear  on  the  bottom  of  the  slab  and  extend  upwards  to  the 
neutral  axis  of  the  slab.  In  the  case  of  joints,  4SIGHT  assumes  that  the  joint  extends  through 
the  entire  depth  of  the  roof  slab.  The  joint  is  filled  with  a joint  compound  with  a known 
permeability  and  service  life. 

The  permeability  of  a cracked  slab  is  calculated  assuming  the  crack  walls  are  smooth 
and  parallel.  Given  a square  slab  with  width  L and  depth  D having  cracks  with  width  w 
penetrating  the  full  depth  D of  the  slab,  the  permeability  of  the  slab  is  a weighted  sum  of 
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(20) 


the  permeability  of  the  crack,  kc , and  the  permeability  of  the  uncracked  concrete,  kQ: 

k c 


w j L 

kc  + ~ — : — k0 


L + w " ' L + w 
•2  L 


w w 


+ 


L -f  w 12  L T w 


■k0 


(21) 


Since  w3  / 12  is  typically  far  greater  than  Lka , the  permeability  of  the  slab  can  be  approxi- 
mated by  the  permeability  due  to  the  crack.  Further,  if  each  of  the  cracks  of  width  w are 
spaced  a distance  a apart,  the  permeability  of  the  slab,  ks,  is 


ks 


(22) 


Joints  can  be  handled  in  a similar  manner  as  cracks.  However,  joints  will  typically  be 
very  much  wider  than  cracks  Since  joints  will  presumably  extend  the  entire  thickness  of  the 
slab,  once  the  joint  fails,  the  flow  through  the  joint  would  overwhelm  the  transport  of  ions 
through  the  central  portion  of  the  slab.  In  fact,  the  transport  coefficients  could  be  as  great 
as,  or  greater  than,  those  of  the  soil.  Therefore,  upon  failure  of  the  joint,  4SIGHT  assumes 
that  the  roof  fails  to  impede  the  flow  of  water  into  the  vault,  the  transport  properties  of 
the  concrete  should  be  approximated  by  the  transport  properties  of  soil,  and  the  calculation 
ceases. 


8 Chemical  Equilibrium 

After  each  time  step,  each  computational  element  is  brought  to  chemical  equilibrium  by 
satisfying  two  conditions: 

1.  If  a salt  exists  as  solid,  the  constituent  ion  concentration  product  equals  the  solubility 
product. 

2.  The  sum  of  the  free  charges  from  all  available  ions  equals  zero,  insuring  local  charge 
neutrality. 

Given  the  salt  ApCa,  composed  of  anion  Aa~  and  cation  C/3+,  condition  1 above  implies 
that  if  ApCa  exists  as  a solid  then 

= K>P  (23) 

where  K3p  is  the  ion  solubility  product.  Condition  2 implies  that  with  m anions  and  n 
cations  present  in  the  pore  solution: 

ra  n 

= £/?cf'  (24) 

«=i  j=i 
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4SIGHT  satisfies  eqns.  23  and  24  using  an  iterative  process.  Each  iteration  is  composed 
of  two  steps: 

1.  Determine  how  many  moles  of  salt  should  be  precipitated  or  dissolved  to  satisfy  eqn.  23. 

2.  Given  K3p  for  water,  what  concentration  of  hydroxyl  ions  ( 0H~ ) is  needed  to  satisfy 
eqn.  24 

These  two  steps  are  repeated,  as  necessary,  until  both  equations  are  satisfied. 

9 Failure  Criteria 

Once  the  rules  for  propagating  ions  through  the  slab  have  been  defined,  criteria  for  terminat- 
ing the  calculations  are  needed  for  each  degradation  mechanism.  Termination  occurs  when 
suitable  failure  criteria  have  been  met. 

9.1  Reinforcement  Corrosion 

Corrosion  of  the  rebar  is  due  to  the  presence  of  chloride  ions  and  proceeds  in  two  stages: 
initiation  and  corrosion.  The  initiation  stage  is  the  time  during  which  chloride  ions  are 
diffusing  through  the  concrete.  During  this  time  there  is  an  insufficient  concentration  of 
chloride  ions  at  the  rebar  for  corrosion  to  occur.  As  the  concentration  of  chloride  ions  at  the 
rebar  increases,  the  pH  decreases  to  insure  charge  neutrality.  When  a sufficient  concentration 
of  chloride  ions  reaches  the  rebar,  corrosion  begins.  The  corrosion  stage  begins  at  the  onset 
of  corrosion  and  lasts  until  failure  of  the  steel  reinforcement,  which  is  typically  only  a few 
years  hence.  Since  the  initiation  stage  may  last  for  hundreds  of  years,  the  duration  of  the 
corrosion  stage  is  insignificant  to  the  total  lifetime  of  the  reinforcement.  Therefore,  an 
accurate  estimate  of  the  service  life  of  the  concrete  can  be  approximated  from  the  duration 
of  the  initiation  period. 

9.2  Sulfate  Attack 

As  the  sulfate  front  proceeds  into  the  concrete,  the  effective  thickness  of  the  concrete  de- 
creases. At  some  point  in  time  the  roof  has  an  insufficient  thickness  to  support  its  load.  At 
this  time,  structural  failure  of  the  roof  will  occur.  This  critical  depth  of  sulfate  penetra- 
tion can  only  be  determined  through  detailed  structural  analysis.  Therefore,  it  is  the  user’s 
responsibility  to  provide  this  information. 
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9.3  Joint  Failure 


The  failure  criteria  for  joints  is  based  upon  the  service  life  of  the  joint  compound.  Since  joint 
widths  will  be  centimeters  wide,  as  previously  stated,  upon  joint  failure  the  permeability  of 
the  roof  slab  will  be  overwhelmed  by  the  flow  through  the  joints.  Therefore,  the  calculations 
terminate  upon  when  the  internal  time  variable  has  reached  the  limit  of  the  joint  service  life. 


10  Calculated  Material  Properties 

All  of  the  concrete  physical  parameters  ( e.g .,  diffusivity,  permeability,  etc.)  are  user  specified 
inputs  to  4SIGHT.  However,  in  cases  where  not  all  properties  are  available,  missing  quantities 
must  be  approximated  using  existing  correlations.  The  physical  properties  must  be  estab- 
lished due  to  both  hydration  and  leaching.  The  physical  properties  due  to  hydration  are  the 
initial  conditions.  However,  as  the  porosity  changes  due  to  leaching,  corrected  values  of  the 
physical  parameters  are  needed. 


10.1  Hydration 


The  hydration  of  cement  can  be  approximated  by  a reaction  between  tri-calcium  silicate 
(C3S)  and  water,  forming  a calcium  silicate  hydrate  (CSH).  A more  elaborate  model  incor- 
porating multiple  mineral  phases  would  require  chemical  analysis  of  the  cement  and  yield 
relatively  little  additional  information  concerning  degree  of  hydration.  The  weight  ratio  of 
water  to  cement,  — , is  the  oft  reported  quantity  to  characterize  the  concrete  mix.  After 
some  period  of  hydration,  the  fraction  of  the  initial  C3S  which  has  hydrated  is  the  degree 
of  hydration,  a.  The  relation  between  these  two  properties  and  porosity  can  be  determined 
stoichiometrically  [7], 


1 + 1.31a 
1 + 3.2- 

C 


(25) 


and  is  valid  for  ^ values  used  in  practice. 

The  diffusivity  can  be  related  to  either  j or  <j>.  After  the  first  100  days  of  hydration  the 
transport  properties  of  most  cement  pastes  are  near  their  asymptotic  values.  Although  the 
values  are  still  changing  after  100  days,  these  changes  are  small  compared  to  the  accuracy 
with  which  these  transport  measurements  can  be  made.  Due  to  this  apparent  steady  state, 
an  empirical  relation  between  Do,  the  chloride  diffusivity  in  m2/s,  and  f was  developed  for 
cement  paste  by  Atkinson,  Nickerson,  and  Valentine[8]  and  Walton,  Plansky,  and  Smith[9] 


w 

log10  D0  = 6.0 9.84 

c 


(26) 


for  values  of  j in  the  range  (0.2-0. 6). 
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D0  for  cement  paste  can  also  be  related  to  (f>  using  the  universal  relation  for  the  formation 
factor  due  to  hydration,  d,  [10]: 

d = = 0.001  + .07(f)2  + 1.8 {(f)  - A8)2H{(f)  - .18)  (27) 

Dci- 

where  D0  is  the  concrete  chloride  diffusivity,  Dj  is  the  free  ion  diffusivity  of  chloride  ions, 
and  H(x)  is  the  Heaviside  function.  The  quantity  7)  is  a constant  for  all  ions. 

The  above  equations  relate  the  diffusivity  of  cement  paste  to  water:cement  ratio  or  poros- 
ity. A relationship  is  now  needed  between  cement  paste  diffusivity  and  concrete  diffusivity. 
Experimental  results  of  Luping  and  Nilssonfll]  for  cement  paste  and  mortar  suggest  that 
their  diffusivities  are  approximately  equal.  Additionally,  results  from  numerical  experiments 
by  Garboczi,  Schwartz,  and  Bentz[12]  investigating  the  effects  of  aggregate-paste  interfacial 
zone  diffusivity  upon  bulk  diffusivity  suggest  that  concrete  diffusivity  is  approximately  equal 
to  the  paste  diffusivity. 

Given  the  permeability  can  be  approximated  from  the  data  in  Hearn,  et  al. [13] : 

k = lO5'0^  x lO-21  m2  (28) 

for  j in  the  range  (0.35,0.80). 

10.2  Leaching 

Once  D and  k have  been  established,  changes  due  to  leaching  can  be  calculated  from  d((f)0), 
where  (f)Q  is  the  initial  porosity  due  to  hydration.  As  the  Ca(OH) 2 is  leached,  the  porosity 
increases.  Let  the  value  of  porosity  after  leaching  be  (f>' , and  the  diffusivity  be  D' . Unfor- 
tunately, the  ratio  D' / D?  is  not  simply  d{<f)')  because  as  calcium  hydroxide  is  leached  from 
the  paste  the  ratio  D' / D*  does  not  retrace  eqn.  27.  Rather,  the  ratio  D' / D*  is  greater  than 
d((f)'),  as  demonstrated  by  the  NIST  microstructural  model. 

The  NIST  cement  microstructural  model[14]  was  used  to  determine  the  relation  for  D' / D0 
upon  leaching.  Results  for  a -^=0.35  paste  are  shown  in  Figure  2.  The  formation  factor 
decreases  with  decreasing  porosity  due  to  hydration,  denoted  by  circles.  Upon  leaching  of 
the  calcium  hydroxide,  the  formation  factor  follows  the  curve  denoted  by  squares.  Given  the 
following  definitions: 

do  = d(<f>0)  d'  = d{<ff)  (29) 

An  empirical  relation  was  developed  to  relate  the  leached  pore  structure  to  the  undamaged 
pore  structure. 

£ = tfo  + 2.0(0' - do)  (30) 

and  is  shown  by  the  solid  curve  in  Figure  2.  Therefore,  the  ratio  of  the  leached  value  of 
diffusivity,  D' , to  the  initial  value  D0  is 

K = 1 

D0  'do 
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Figure  2:  Formation  factor,  D / D* , due  to  hydration  (circles)  and  to  leaching  (squares)  for 
w/c=0.35  (after  [14]).  The  solid  line  represents  the  approximation  used  by  4SIGHT  for  the 
leaching  formation  factor. 

The  relative  change  in  permeability  can  be  calculated  using  £.  The  Katz-Thompson 
equation  relates  permeability  to  d [15]: 


d 2 

k = -±d 
226 


(32) 


where  dc  is  the  diameter  of  the  largest  sphere  which  can  pass  through  the  pore  space  of  the 
sample.  Also,  given  that  dc  oc  d [16],  the  relative  change  in  permeability  becomes 


kQ  l 


(33) 


The  final  dimensionless  advection-diffusion  equation  is  a combination  of  eqns.  14,  31,  and 


33: 


S=v-lvp+(!)  vd'Vc- 


(34) 


There  is  one  such  equation  for  chemical  species  i. 


11  Program  Outline 

The  flow  of  the  program  is  summarized  by  the  pseudo-code  program  shown  in  Figure  3.  The 
program  naturally  divides  itself  into  three  parts:  initialization,  ion  propagation,  and  final 
state  output. 
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Figure  3:  4SIGHT  flow  chart. 
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12  Assumptions 

The  calculations  of  4SIGHT  are  based  upon  a set  of  assumptions.  For  completeness,  these 
are  assumptions  are  enumerated  here: 

*v 

1.  The  vault  system  is,  and  remains,  water  saturated.  This  should  be  valid  for  all  but  the 
most  arid  locations. 

2.  Darcy  flow  is  valid,  even  at  the  bottom  of  the  slab  which  might  not  be  in  direct  contact 
with  liquid. 

3.  There  is  sufficient  oxygen  present  for  corrosion.  This  a conservative  approach  to  cor- 
rosion. 

4.  The  concrete  vault  begins  service  approximately  100  days  after  casting.  Therefore,  the 
transport  properties  have  reached  about  90%  of  their  assymptotic  values. 

5.  The  external  ion  concentration  and  hydraulic  head  are  constant  over  the  life  of  the 
vault. 

6.  As  the  sulfate  front  advances,  the  concrete  in  its  wake  can  be  treated  like  soil  and  the 
external  conditions  can  be  advanced  to  this  point. 
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Appendix  A 


A Input  Parsing 

The  input  parsing  was  accomplished  using  public  domain  versions  of  LEX  and  YACC  written 
for  the  PC.  These  tools  simplify  parsing  user  input,  facilitating  modification  and  extensions 
in  program  capabilities.  The  user  unfamiliar  with  LEX  and  YACC  should  read  the  book  by 
Mason  and  Brown  listed  in  the  references.  To  obtain  a copy  of  the  public  domain  version, 
access  the  following  directory  through  anonymous  ftp: 

machine:  omnigate.clarkson.edu 
directory:  pub/msdos/djgpp 

LEX  and  YACC  are  programs  that  take  as  input  user  specifications  for  a particular  syntax 
and  create  a C file  which  performs  the  functions  specified  in  the  user  LEX  and  YACC  files. 
The  resulting  C file  is  then  simply  #include-ed  into  the  main  C program.  Of  course,  the 
user  can  code  the  parsing  directly  in  C,  but  LEX  and  YACC  make  the  coding  less  errorsome 
and  more  extensible. 


A.l 


A.l  LEX  Specification 


•/.{ 

#include 

"4sight_t .h" 

•/.> 

integer 

[0-9]  + 

dreal 

([0-9]+"."  [0-9]*)  | ([0-9]*"  [0-9]+) 

ereal 

({dreal} I {integer}) [eE] [+-] ? [0-9] + 

real 

{dreal} I {ereal} 

diff 

[Dd]  [Ii]  [Ff]  [Ff] 

perm 

[Pp]  [Ee]  [Rr]  [Mm] 

wc 

[Ww]  [Cc] 

thick 

[Tt]  [Hh]  [Ii]  [Cc]  [Kk]  [Nn]  [Ee]  [Ss]  [Ss] 

neuax 

[Nn]  [Aa]  [Xx]  [Ii]  [Ss] 

youngs 

[Yy]  [Oo]  [Uu]  [Nn]  [Gg]  [Ss] 

beta 

[Bb]  [Ee]  [Tt]  [Aa] 

ce 

[Cc]  [Ee] 

rough 

[Rr]  [Oo]  [Uu]  [Gg]  [Hh]  [Nn]  [Ee]  [Ss]  [Ss] 

gamma 

[Gg]  [Aa]  [Mm]  [Mm]  [Aa] 

poisson 

[Pp]  [0o]  [Ii]  [Ss]  [Ss]  [Oo]  [Nn] 

extrn 

[Ee]  [Xx]  [Tt]  [Ee]  [Rr]  [Nn]  [Aa]  [LI] 

intrn 

[Ii]  [Nn]  [Tt]  [Ee]  [Rr]  [Nn]  [Aa]  [LI] 

time 

[Tt]  [Ii]  [Mm]  [Ee] 

depth 

[Dd]  [Ee]  [Pp]  [Tt]  [Hh] 

rebar 

[Rr]  [Ee]  [Bb]  [Aa]  [Rr] 

output 

[Oo]  [Uu]  [Tt]  [Pp]  [Uu]  [Tt] 

head 

[Hh]  [Ee]  [Aa]  [Dd] 

crack 

[Cc]  [Rr]  [Aa]  [Cc]  [Kk] 

j oint 

[Jj]  [0o]  [Ii]  [Nn]  [Tt] 

at 

[Aa]  [Tt] 

until 

[Uu]  [Nn]  [Tt]  [Ii]  [LI] 

assign 

ii _ ii  | ii  . ii 

wspace 

[ \t]  + 

kill 

( ( [Qq]  ( [Uu]  [Ii]  [Tt]  ) ?)  | ( [Ee]  ( [Xx]  [Ii]  [Tt]  ) ?)  ) 

nl 

\n 

•///. 

{wspace} 


OH 

{ 

return 

HYDROXIDE;} 

Cl 

{ 

return 

CHLORINE;} 

S04 

{ 

return 

SULFATE;} 

A. 2 


C03 

{ 

return  CARBONATE;} 

H 

{ 

return  HYDROGEN;} 

Ca 

{ 

return  CALCIUM;} 

Na 

{ 

return  SODIUM;} 

K 

{ 

return  POTASSIUM;} 

{diff} 

{ 

return  DIFF;} 

{perm} 

{ 

return  PERM;} 

{wc} 

{ 

return  WC;} 

{thick} 

{ 

return  THICKNESS;} 

{neuax} 

{ 

return  NEUAX;} 

{youngs} 

{ 

return  YOUNGS;} 

{beta} 

{ 

return  BETA;} 

{ce} 

{ 

return  Ce;} 

{rough} 

{ 

return  ROUGHNESS;} 

{gamma} 

{ 

return  GAMMA;} 

{poisson} 

{ 

return  POISSON;} 

{extrn} 

{ 

return  EXTRN ; } 

{intrn} 

{ 

return  INTRN;} 

{time} 

{ 

return  TIME;} 

{depth} 

{ 

return  DEPTH;} 

{rebar} 

{ 

return  REBAR;} 

{output} 

{ 

return  OUTPUT;} 

{head} 

{ 

return  HEAD;} 

{crack} 

{ 

return  CRACK;} 

{joint} 

{ 

return  JOINT;} 

{at} 

{ 

return  AT;} 

{until} 

{ 

return  UNTIL;} 

{integer} 

{ 

sscanf  (yytext , '"/.Id" , &yylval . int32)  ; 

return  INT32;} 

{real} 

{ 

sscanf  (yytext , "'/.If" , &yylval . real)  ; 

return  REAL;} 

{assign} 

{ 

return  EQUALS;} 

{nl} 

{ 

return  NEWLINE;} 

{kill} 

{ 

return  KILL;} 

• 

{ 

return  yytext[0];} 

XX 

A. 3 
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A. 2 YACC  Specification 
U 

int  yylex(void) ; 
void  yyerror(char  *)  ; 
void  yyerror(s) 
char  *s; 

{ 

/*  fprintf  (stderr , "'/,s  \t  see  line  '/.d\n"  , s , input_f  ile.line) ; */ 

> 

void  syntax_error (char  *) ; 
void  syntax_error (char  *s) 

fprintf  (stderr , "'/,s  \t  see  line  '/.d\n"  , s , input _file_line)  ; 

> 

’/.} 

'/.union  { 

double 
long 
int 

> 

'/.token  <real>  REAL 
'/.token  <int32>  INT32 
'/.token  KILL 
'/.token  EQUALS 
'/.token  NEWLINE 
'/.token  DIFF 
'/.token  PERM 
'/.token  WC 
'/.token  THICKNESS 
'/.token  NEUAX 
'/.token  YOUNGS 
'/.token  BETA 
'/.token  Ce 
'/.token  ROUGHNESS 
‘/.token  GAMMA 


real;  /*  real  value  */ 

int32;  /*  DoubleWord  Integer  */ 

int 16;  /*  SingleWord  Integer  */ 


A. 4 


’/.token  POISSON 
•/.token  EXTRN 
•/.token  INTRN 
•/.token  HYDROGEN 
•/.token  CALCIUM 
•/.token  SODIUM 
‘/.token  POTASSIUM 
•/.token  HYDROXIDE 
‘/.token  CHLORINE 
•/.token  SULFATE 
‘/.token  CARBONATE 
•/.token  TIME 
‘/.token  DEPTH 
'/.token  REBAR 
‘/.token  OUTPUT 
'/.token  HEAD 
'/.token  CRACK 
‘/.token  JOINT 
'/.token  AT 
‘/.token  UNTIL 


'/.type  <real>  expr 
V.type  <intl6>  anion 
'/.type  <int!6>  cation 


•/.left  • + > 
'/.left  '*»  >/> 

*/.*/. 


lines : 


/*  NOTHING  */ 

lines  line 


line:  NEWLINE  {input_f ile_line++ ; } 

I terminate  { printf("'/,d  lines  parsed . \n"  , input _file_line) 
return  0;} 

I assign 
I recvr  NEWLINE 
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terminate 


recvr : 


expr : 


assign : 


KILL  {printf ("Received  KILL.\n");} 


error  { syntax.error ("Syntax  error");} 

recvr  error 

recvr  EQUALS  expr  { syntax_error ("Unknown  Command");} 


INT32  {$$  = (real) $1 ; } 

REAL  {$$  = $1;} 


DIFF  EQUALS  expr  { 
PERM  EQUALS  expr  { 
WC  EQUALS  expr  { 
THICKNESS  EQUALS  expr  { 

NEUAX  EQUALS  expr  { 
YOUNGS  EQUALS  expr  { 
BETA  EQUALS  expr  { 
Ce  EQUALS  expr  { 
ROUGHNESS  EQUALS  expr  { 
GAMMA  EQUALS  expr  { 
POISSON  EQUALS  expr  { 


EXTRN  cation  EQUALS  expr 
EXTRN  anion  EQUALS  expr 
INTRN  cation  EQUALS  expr 


Dinfty. value  = $3; 

Dinf ty . is.default  = FALSE;} 

kinfty. value  = $3; 

kinfty . is_default  = FALSE;} 

wc. value  = $3; 

wc . is.def ault  = FALSE;} 

thickness .value  = $3; 

thickness . is_default  = FALSE; 

sample.length  = thickness . value ; } 

neutral_axis_depth . value  = $3; 

neutral_axis_depth . is_def ault=FALSE; } 

Youngs. value  = $3*1.0E+09; 

Youngs . is_default  = FALSE;} 

beta. value  = $3; 

beta. is_def ault  = FALSE;} 

CE. value  = $3; 

CE . is_def ault  = FALSE;} 

roughness .value  = $3; 

roughness . is_def ault  = FALSE;} 

gamma. value  = $3; 

gamma. is .default  = FALSE;} 

nu. value  = $3; 

nu . is.def ault  = FALSE;} 

{ cation [$2] . c [0]  = $4;} 

{ anion  [$2]  .c[0]  = $4;} 

{ cation [$2] . c [1]  = $4;} 
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DEPTH  EQUALS  expr 

REBAR  EQUALS  expr 

OUTPUT  cation  anion 
OUTPUT  cation 
OUTPUT  anion 
HEAD  EQUALS  expr 


I INTRN  anion  EQUALS  expr  { anion[$2] .c[l]  = $4;} 

I TIME  EQUALS  expr  { terminat ion.type  = TIME.LIMIT; 

if (MaxDay . is_def ault==TRUE) 

MaxDay. value  = $3; 
else 

MaxDay. value  = MIN (MaxDay .value ,$3) ; 
MaxDay . is_default  = FALSE;} 

{ termination_type  = STRUCT_LIMIT; 

MaxDepth  = $3;} 

{ rebar_depth=$3 ; 

terminat ion_type  = STRUCT_LIMIT; } 

{ sol_array [$2] [$3] .output_flag=TRUE;} 

{ cat ion [$2] . output _f lag  = TRUE;} 

{ anion [$2] .output_flag  = TRUE;} 

{ head. value  = $3; 

head . is_def ault  = FALSE;} 

CRACK  EQUALS  expr  AT  expr  DEPTH  expr 

{ crack.width. value  = $3; 

crack.width. is_def ault  = FALSE; 
crack_spacing  = $5; 
crack_depth  = $7;} 

JOINT  EQUALS  expr  AT  expr  UNTIL  expr 
{ joint.width  = $3; 
j oint_spacing  = $5; 
joint_lif etime  = $7*365; 
if (MaxDay . is_def ault  ==  TRUE) 

MaxDay. value  = j oint_lif etime ; 
else 

MaxDay. value  = MIN(MaxDay. value, j oint_lif etime) ; 
MaxDay . is_def ault  = FALSE; 
joint_is_specif ied  = TRUE;} 

JOINT  PERM  EQUALS  expr  { j oint.permeability  = $4;} 


cation:  HYDROGEN 

I CALCIUM 
I SODIUM 
I POTASSIUM 


{ $$ 
{ $$ 
{ $$ 
{ $$ 


H;} 

Ca;} 

Na;} 

K;} 
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anion: 

HYDROXIDE 

{ 

$$  = OH;} 

1 CHLORINE 

{ 

$$  = Cl;} 

I SULFATE 

{ 

$$  = S04 ; } 

I CARBONATE 

{ 

$$  = C03 ; } 
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B CWEB  - Source  Code 


The  source  code  (except  for  the  parsing  described  above)  was  written  using  the  literate 
programming  tool  CWEB  written  by  Donald  Knuth  and  Silvio  Levy.  The  CWEB  tool  is  composed 
of  two  programs:  cweave  and  ctangle.  The  user  writes  a CWEB  file  using. any  ASCII  editor. 
This  file  contains  both  the  C code  and  the  documentation,  ctangle  takes  this  input  file 
and  extracts  the  C code.  Similarly,  cweave  takes  the  input  file  and  creates  a T^Xfile  which 
can  be  Tj^X-ed,  creating  a formatted  version  of  the  complete  documentation,  including  the 
source  code.  Therefore,  to  use  these  tools,  the  user  must  have  both  CWEB  and  T^Xfor  the 
PC.  Fortunately,  both  of  these  tools  are  public  domain,  and  can  be  accessed  by  anonymous 
ftp: 


TeX 

machine: 

ftp . nj it . edu 

directory: 

pub/msdos/ emtex 

CWEB 

machine: 

labrea. Stanford. edu 

directory: 

pub/cweb 

The  only  modification  made  to  CWEB  for  the  development  of  4SIGHT  was  the  extension 
to  LTgX.  To  use  the  extensions  of  LTgXwith  CWEB,  one  only  needs  a LTgXstyle  sheet.  This 
too  is  available  by  anonymous  ftp: 

machine:  ftp.th-darmstadt.de 

directory:  pub/programming/literate-programming/c . C++ 

The  pages  that  follow  are  verbatim  output  from  cweave,  via  LTgX: 
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1.  INTRODUCTION. 

The  4SIGHT  computer  program  has  been  developed  by  the  Inorganic  Group  of  the  Building  Mate- 
rials Division  at  NIST.  This  program  facilitates  estimating  the  service  life  of  low  level  nuclear  waste 
storage  sites  by  determining  the  hydraulic  conductivity  of  the  vault  roof  as  a function  of  time.  This 
objective  is  acheived  by  considering  prominent  deterioration  mechanisms  to  predict  the  intrinsic 
permeability,  effluent,  and  working  thickness  of  the  vault  slab. 

The  condition  of  the  vault  slab  is  assesssed  by  modelling  the  ingress  of  deleterious  ion  species 
by  both  diffusion  and  convective  flow  due  to  hydraulic  pressure  gradients.  Additionally,  changes  in 
physical  properties  due  to  the  leaching/deposition  of  salts  are  incorporated  into  the  program.  The 
propagation  of  ions  are  calculated  using  an  advection-diffusion  equation.  After  each  discrete  time 
step,  each  computational  cell  is  brought  to  chemical  equilibrium  through  precipitation/dissolution 
of  ionic  salts. 

The  program  incorporates  a one-dimensional  model  of  the  vault  since  the  critical  deterioration 
will  most  likely  occur  at  the  ceiling  of  the  vault.  The  ceiling  is  the  most  likely  element  to  experience 
moisture,  stress,  and  external  ions  such  as  Cl~  and  SO4.  Addtionally,  the  majority  of  the  transport 
will  occur  at  the  center  of  the  roof  slab  since  the  edges  of  the  slab  will  be  supported.  Therefore,  a 
simple  one-dimension  model  of  the  flow  at  the  center  of  the  slab  will  give  a conservative  estimate  of 
the  transport  for  the  entire  slab. 

The  program  proceeds  in  a strainghtforward  manner;  The  program  gets  input  parameters  from 
the  user;  initializes  the  ion  information  and  scaling  parameters;  iterates  diffusion  and  convection 
calculations,  bringing  each  computational  element  to  chemical  equilibrium  between  iterations;  de- 
termines whether  failure  has  occured;  and  prints  pertenent  information  about  the  performance  of 
the  vault. 

The  output  of  the  program  can  be  printed  to  either  the  screen  or  to  an  ASCII  file.  The  ASCII 
file  can  then  be  used  as  input  to  a spreadsheet  program  for  further  analysis. 


' 
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2.  MAIN. 

The  body  of  4SIGHT  is  fairly  straightforward  and  the  logical  flow  can  be  inferred  from  the  outline 
below. 

#define  banner  "\"\tThisuisu4SIGHTu(Versionul . 0)\"\n\n" 

#define  NUM.CELLS  20 
# define  NUM.SURFACES  (NUM.CELLS + 1) 

# define  DLIMIT  0.0005 
( Include  files  3 ); 

(Preprocessor  definitions  li); 

(User  defined  data  types  14); 

(Global  variables  4); 

(Function  declarations  23); 

(Interrupt  handlers  24); 

(Input  parsing  routine  13); 
int  main( int  argc,  char  **argv) 

{ 

int  i,  j,  k; 

( System  setup  5 ); 

( Parse  input  6 ); 

(Parameter  initialization  7); 

( Print  header  75 ); 

(Print  intermediate  results  76); 
do  { 

if  ( kbhit( ))  { 
i = getch( ); 
control-break ( ); 

} 

( Ion  transport  8); 

(Adjust  parameters  9); 

(Print  intermediate  results  76); 

(Assess  termination  10); 

} while  (-> termination)-, 

(Print  termination  information  77); 

(Print  desired  output  78); 
return  1; 

} 


3.  The  following  are  #include  files  needed  for  some  of  the  intrinsic  functions  such  as  pnntf  and 
scanf . The  math  functions  pow  and  log  and  the  math  error  routine  matherr  require  math.h. 
Additional  include  files  are  declared  when  needed.  The  routine  kbhit()  requires  the  header  file 
conio.h. 

( Include  files  3)  = 

^include  <stdio.h> 

#include  <math.h> 

^include  <conio.h> 

^include  <bios.h> 

^include  <time.h> 

See  also  sections  22,  28,  and  32. 

This  code  is  used  in  section  2. 
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4.  ( Global  variables  4 ) = 

long  elapsed-time ; 

See  also  sections  12,  19,  20,  21,  29,  31,  40,  63,  and  69. 
This  code  is  used  in  section  2. 


5.  System  Setup. 

#define  CTRL_BRK_0N  1 
#define  GET.TIME  0 
( System  setup  5 ) = 
print/  (banner)] 

elapsedJime  = &iositme(GET_TIME, 0 l);  /*  get  BIOS  timer  value  */ 

(Determine  STDIN  30); 

i = se/c6r£(CTRL_BRK_0N);  /*  check  ctrl-brk  every  system  call  */ 
ctrlbrk  (control-break); 
signal( SIGFPE,  div.O ); 
initialize-ion.data  ( ); 

This  code  is  used  in  section  2. 


6.  Parse  Input. 

( Parse  input  6 ) = 
i = yyparse( ); 

This  code  is  used  in  section  2. 


7.  Parameter  Initialization. 

(Parameter  initialization  7)  = 
initialize-parameters( ); 
update-pressures ( ); 
aci-guidlines( ); 
chemicaLequilibrium(F  ALSE); 
cheTnical-equilibrium(E  ALSE); 
chemical-equilibrium  (FALSE); 

( Print  initial  system  state  74 ); 

This  code  is  used  in  section  2. 


8.  Ion  Transport. 

( Ion  transport  8 ) = 
ion-diffusion ( ); 
chemicaLequihbnum(  TRUE); 
advancesulfaie-front( ); 

This  code  is  used  in  section  2. 


9.  Adjust  Parameters. 

(Adjust  parameters  9)  = 

(Advance  global  clock  70 ); 
(Adjust  physical  parameters  68); 
update-pressures  ( ); 

(Update  the  porosities  67); 

This  code  is  used  in  section  2. 
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10.  Assess  Termination. 

(Assess  termination  10)  = 

(Assess  chloride  penetration  64); 
(Assess  simulation  termination  72); 
(Check  time  dependencies  7i); 

This  code  is  used  in  section  2. 


11.  The  following  are  useful  macro  definitions  which  save  a lot  of  typing. 

(Preprocessor  definitions  ll ) = 

#define  SQR(a)  ((a)  * (a)) 

#define  CUB(a)  ((a)  * (a)  * (a)) 

#define  MAX(a,6)  ((a  > 6)  ? a : b) 

#define  HIN(a,6)  ((a  < b)  ? a : b) 

See  also  section  46. 

This  code  is  used  in  section  2. 


12.  Declare  yyparse  which  is  written  in  YACC  and  LEX.  The  function  yyparse  takes  no  arguments 
and  simply  parses  the  input,  assigning  values  to  global  variables.  The  routine  does  not  “return” 
until  the  entire  input  stream  has  been  parsed. 

(Global  variables  4)  -f= 

int  input.fileJine  = 1;  /*  line  number  of  input  file  (lex/yacc)  */ 


13.  (Input  parsing  routine  13)  = 
int  yj/parse(void); 

^include  "yyparse. c" 

This  code  is  used  in  section  2. 
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14.  User  defined  data  structures. 

These  are  constructed  datatypes  used  by  4SIGHT.  Two  fundamental  types  are  real  and  boolean. 
Using  real  lets  the  programmer  easily  change  between  float  and  double  data  types  for  real  numbers. 

The  boolean  types  are  defined  {FALSE, TRUE}  because  Borland  assigns  the  first  entry  the  number 
0.  Variables  defined  as  boolean  can  be  used  in  a conditional  statement  as  i{(myboolvar ) without 
direct  comparison  to  FALSE  or  TRUE. 

(User  defined  data  types  14)  = 
typedef  double  real; 
typedef  enum  { 

FALSE, TRUE 
} boolean; 

See  also  sections  15,  16,  and  17. 

This  code  is  used  in  section  2. 


15.  The  ION  data  type  contains  all  important  information  about  a single  ion.  This  information 
includes: 

c[j:  concentration  at  each  computational  cell  surface. 
tj:  the  ratio  Dji/ D0  for  each  ion. 

Dj:  free  ion  diffusivity. 
valence : the  valence  of  the  ion. 

output-flag : a flag  to  denote  whether  the  final  ion  concentration,  as  a function  of  distance,  is  printed 
in  system. out. 

name []:  the  ASCII  name  of  the  ion. 

(User  defined  data  types  14)  += 
typedef  struct  { 

real  c[NUM_SURFACES];  /*  concentration  at  each  surface  */ 

real  moles  [NUM_SURFACES];  /*  moles  ion  in  solution  */ 

real  77;  /*  Dfi/Do  */ 

real  Dj\  /*  ion  free  diffusivity  */ 

real  valence ; /*  ion  valence  */ 

boolean  output-flag ; /*  print  to  system. out?  */ 

char  name  [5];  /*  ASCII  name  of  ion  */ 

} ION; 

16.  The  SOLARRAY  data  type  contains  the  solubility  products  and  stiochiometric  ratios  for  each 
anion-cation  combination.  The  salt  is  CmAn,  where  C represents  the  cation,  and  A the  anion.  The 
quantity  of  solid  is  defined  on  the  cell  surfaces. 

(User  defined  data  types  14)  += 
typedef  struct  { 

real  s[NUM_SURFACESj;  /*  moles  solid  salt  in  element  */ 

real  ksp;  /*  salt  solubility  constant  */ 

real  m;  /*  stiochiometric  cation  factor  */ 

real  n;  /*  stiochiometric  anion  factor  */ 

real  molar-density ; /*  cm2 /mole  */ 

boolean  output.flaq:  /*  print  salt  to  system. out?  */ 

} SOLARRAY; 
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17.  The  MATERIAL  datatype  is  used  for  any  material  property  that  can  be  specified  by  the  user. 
It  carries  a flag  to  denote  whether  the  value  is  the  default  one,  or  the  one  specified  by  the  user. 


(User  defined  data  types  14)  += 
typedef  struct  { 
real  value ; 
boolean  is.default ; 

} MATERIAL; 
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18.  Global  variables. 

Most  of  these  are  variables  that  will  be  changed  by  the  input  file.  The  listing  is  broken  into  major 
groups:  material  parameters,  scaling  parameters,  etc. 


19.  Transport  parameters.  These  are  the  parameters  defined  at  each  element  that  determine  the 
rate  of  transport.  These  include  the  inverse  formation  factor  £,  pressure  't,  and  the  porosity  </>. 
(Global  variables  4)  += 

real  £[NUM_CELLS];  /*  changing  formation  factor  */ 

real  £'[NUM_SURFACES];  /*  changing  formation  factor  at  the  surface  */ 

real  ^[NUH.SURFACES];  /*  Pressure  */ 

real  </>n(NUM_CELLS];  /*  current  porosity  */ 

real  </>,n[NUM_SURFACES];  /*  porosity  at  the  surface  */ 

real  (/>/n-1[NUM_SURFACES];  /*  previous  porosity  at  concentration  */ 

real  /ilre[NUM_SURFACES];  /*  litres  of  solution  (surface  area)*/ 

real  strength  [NUM_SURFACES];  /*  strength  in  psi  */ 

int  FIRST_CELL  = 0;  /*  first  active  cell  number  */ 


20.  User  specified  parameters.  These  parameters  determine  the  ultimate  properties  of  the  sample. 
Although  default  values  are  given,  user  specified  values  will  override  these  values. 

( Global  variables  4 ) += 

MATERIAL  thickness  = {1.0,  TRUE};  /*  roof  thickness  */ 

MATERIAL  neutral-axis. depth  = {0.75,  TRUE};  /*  depth  of  neutral  axis  */ 

MATERIAL  k^  = {2.5-  10-18,TRUE};  /*  ultimate  permeability  */ 

MATERIAL  Doo  = {5.7  ■ 10-12,TRUE};  /*  ultimate  diffusivity  */ 

MATERIAL  head  = {5.0, TRUE};  /*  external  pressure  head  */ 

MATERIAL  wc  = {0.45, TRUE};  /*  water:cement  */ 

MATERIAL  Youngs  = {20.0  • 10+09, TRUE};  /*  Young’s  modulus  */ 

MATERIAL  /3  = {1.8-  10-O6,TRUE};  /*  linear  strain  coefficient  */ 

MATERIAL  CE  = {350., TRUE};  /*  cone,  sulfate  as  ettringite  */ 

MATERIAL  roughness  = {1.0, TRUE};  /*  roughness  factor  */ 

MATERIAL  7 = {10.0, TRUE};  /*  fracture  surface  energy  */ 

MATERIAL  u = {0.2, TRUE};  /*  Poisson  ration  */ 


21.  Scaling  parameters. 

(Global  variables  4)  += 

real  L = 0.010;  /*  length  of  ELEMENTS  so  that  AX  = 1 */ 

real  D0  = 1.0  • 10— 1 1 ; /*  scaling  diffusivity  */ 

real  k0  = 1.0  • 10-18;  /*  scaling  permeability  */ 

real  AT  =1.0;  /*  dimensionless  time  increment  */ 

real  AA[NUM_SURFACES];  /*  dimensionless  length  between  each  cell  */ 
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22.  Interrupt  handlers. 

The  following  routines  handle  ctrl-brk  and  math  exceptions.  The  routine  ctrlbrk  requires  dos.h. 
The  routine  signal  requires  signal. h. 

( Include  files  3)  += 

^include  <dos.h> 

^include  <signal.h> 

23.  (Function  declarations  23)  = 
int  control-break  (\ oid); 

int  matherr (struct  exception  *); 
void  div-0(\ oid); 

See  also  sections  33,  39,  43,  51,  54,  56,  59,  60,  61,  and  65. 

This  code  is  used  in  section  2. 


24.  control-break  handles  interruptions  due  to  the  use  striking  CTRL-BRK  during  execution. 
( Interrupt  handlers  24  ) = 
int  control-break  (void)  { 
int  i,  j,  k ; 

pnntf("\ n\n\tUSERuCTRL-BRK!  .\n\nM);  (Print  desired  output  78) 
return  0;  } 

See  also  sections  25  and  26. 

This  code  is  used  in  section  2. 


25.  The  routine  matherr  handles  exceptions  from  the  math  coprocessor. 

( Interrupt  handlers  24  ) += 

int  matherr  (struct  exception  *e) 

{ 

pnnf/("\  n\n"); 
print/("\tmatherr  :u"); 
switch  (e-type)  { 
case  DOMAIN:  print/(  "DOMAIN"); 
break; 

case  SING:  print/(  "SINGULARITY"); 
break; 

case  OVERFLOW:  print/("0VERFL0W"); 
break; 

case  UNDERFLOW:  print/ ("UNDERFLOW"); 
break; 

case  TLOSS:  print/("TL0SS"); 
break; 

default:  print/ ("UNKN0WNuERR0RuTYPE ! " ); 
break; 

} 

pnntf("\  n"); 
pnntf(‘'\  n\n\t"); 

print/("f  unction:  \t'/,s  ('/, If  ,*/, If  )\n",  e-name , e-aryf  , e-ary2); 

pnntf("\ n"); 

erit(O); 

} 
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26.  The  routine  div.O  handles  possible  divide  by  zero  errors. 

(Interrupt  handlers  24)  += 
void  </iv_0(void) 

{ 

int  i,  j,  k] 

pnn//("\n\n\tPossibleudivideubyuO ! \n\n"); 

(Print  desired  output  78) 
exii(  1); 

} 
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27.  INITIALIZATION. 

28.  Redirection  of  sidin.  Allow  the  user  to  either  use  redirection  at  the  command  line  or  to  specify 
the  input  file  directly.  If  an  input  file  is  given  at  the  command  line  without  redirection  (<),  re-define 
sidin  to  come  from  the  input  file.  The  ^include  files  (f  cntl  .h  and  sys/stat  .h)  are  for  the  routine 
dup2(),  and  errno.h  is  required  for  the  global  variable  errno. 

# define  STDIN  0 
(Include  files  3)  += 

#include  <fcntl.h> 

^include  <sys/stat.h> 

#include  <errno.h> 


29.  If  redirection  is  not  used,  set  input-file-handle  to  the  file  specified  at  the  command  line. 
(Global  variables  4)  += 

int  input-file-handle ; /*  file  containing  input  parameters  */ 


30.  If  an  input  file  was  given  without  redirection,  use  dup2( ) to  copy  the  input  file  handle  number 
to  the  existing  sidin  file  handle. 

(Determine  STDIN  30)  = 
if  ( argc  > 1)  { /*  redirect  from  input  file  */ 

fprintf  (stderr , "Calculating.  . ."); 

if  (( input-file-handle  = open  ( argv  [1],  0_RD0NLY))  = —1)  { 
switch  (errno)  { 

case  ENOENT:  fprintf  (stderr , "Nousuchuf  ile:uy,s\n",  ar<7r[l]); 
break; 

case  EMFILE:  fprintf  (stderr , "Tooumanyuopenuf  iles . \n"); 
break ; 

case  EACCES:  fprintf  (stderr , "Permissionudenied.\n"); 
break; 

case  EINVACC:  fprintf  (stderr , "Invaliduaccessucode.  \n"); 
break ; 

default:  fprintf  (stderr , "Error : uUnknownuerrorijCode . \n"); 
break ; 

} 

exit(  0); 

} 

dup2  (input-file-handle , STDIN); 

} 

else  /*  interactive  mode  */ 
pn7d/("Enteru commands : \n"); 

This  code  is  used  in  section  5. 
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31.  Initialize  solubility  array. 

soLarray  is  vital  to  determining  the  chemical  equilibrium  of  each  element.  This  procedure  also 
initializes  the  two  ION  arrays:  anion  and  cation 
(Global  variables  4)  += 
int  num.cations; 
int  num.anions  \ 

SOLARRAY  **sol.array, 

ION  *cation ; 

ION  *anion ; 


32.  The  routine  malloc()  requires  stdlib.h. 
(Include  files  3)  -f= 

#include  <stdlib.h> 


33.  (Function  declarations  23)  += 
void  initiahze.ion.data  (void); 


34.  Initializing  ion  and  salt  array.  The  data  for  the  ions  and  salts  is  in  the  file  ion.db.  Must  first 
determine  the  number  of  anions  and  cations,  including  OH~  and  H+ , and  then  allocate  memory 
for  anion,  cation,  and  soLarray. 
void  initialize Jon-data(  ) 

{ 

int  i,  j,  k\ 

FILE  *ion.database;  /*  ion  and  salt  data  base  */ 
if  (( ion-database  = fopen(" ion.db",  "rt"))  = A)  { 

print/ ("\ n\t  ERROR : uuUnableutouopenuion . db ! \n\n"); 
eri/(0); 

} 

/scan/ (ion.database , "’/.du'/.d",  Sznum.cations , hnum.anions)-, 

( Allocate  ion  and  salt  arrays  35 ) 

(Input  ion  data  36) 

(Input  salt  data  37) 
fclose(ion.database ); 

(Zero-out  ion  and  salt  data  38) 

} 


35.  Allocate  memory  for  all  of  the  arrays.  The  total  number  of  cations/anions  includes  OH  and 

H+. 

( Allocate  ion  and  salt  arrays  35 ) = 

anion  = (ION  *)  maUoc(num.anions  * sizeof (ION)); 

cation  = (ION  *)  malloc(num.cations  * sizeof  (ION));  /*  allocate  2-D  array  */ 

soLarray  = (SOLARRAY  **)  malloc(num-cations  * sizeof  (SOLARRAY  *)); 
for  (i  = 0;  i < num.cations',  i++)  { 

sol  array[i]  = (SOLARRAY  *)  malloc(num.amons  * sizeof  (SOLARRAY)); 

} 

This  code  is  used  in  section  34. 
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36.  Read  in  the  ion  data  from  ion-database , cations  first,  anions  second. 

(Input  ion  data  36)  = 

for  (i  = 0;  i < num.cahons;  i++)  { 

fscanf  (ion-database  , cation[i].name ); 

fscanf  (ion.database , "'/,lf  ",  &ication[i\.  valence); 
fscanf  (ion-database , "'/,le",  &ccation[i].Dj); 
cation[i]. output-flag  = FALSE; 

} 

for  ( j = 0;  j < num.anions ; j++)  { 

fscanf  (ion-database , anion[j].name); 

fscanf  (ion-database , If  anion  [j]. valence); 
fscanf  (ion-database , "*/,le",  k.anion \j].Dj); 
anion\j], output-flag  = FALSE; 

} 

This  code  is  used  in  section  34. 


37.  Read  salt  data  in  cation  major  order. 

(Input  salt  data  37)  = 

for  (i  = 0;  i < num-cations;  i++) 
for  ( j = 0 ; j < num.anions;  j++)  { 

fscanf  (ion.database , "'/.le",  &soLarray[i][j].ksp); 

fscanf  (ion.database , ",  &tsol.array  [i][j].m); 

fscanf  (ion.database , "'/.If  ",  &csol.array[i][j].n); 

fscanf  (ion.database , "'/.If",  &soLarray[i]\j].molar-density); 

sol-array  [i](j].  output-flag  = FALSE; 

} 

This  code  is  used  in  section  34. 


38.  Initialize  the  ion  and  salt  data  to  zeros. 
( Zero-out  ion  and  salt  data  38 ) = 
for  (Ar  = 0;  k < NUM.SURFACES;  k++)  { 
for  (i  = 0;  i < num.cations;  i++)  { 
cation[i\.c[k\  = 0.0; 
cation[i].moles[k\  = 0.0; 
for  (j  = 0;  j < num.anions;  j++)  { 
anion[j].c[k]  = 0.0; 
anion[j].moles[k]  = 0.0; 
S0/_array[t][j].s[fc]  = 0.0; 

} 

} 

} 

This  code  is  used  in  section  34. 
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39.  Initialize  parameters. 

The  scale  and  physical  parameters  must  be  set 
ion  initialization  routines. 

^define  OH  0 
# define  Cl  1 
^define  S04  2 
#define  C03  3 
# define  H 0 
# define  Ca  1 
# define  Na  2 
# define  K 3 
#define  TIHE.LIMIT  0 
define  STRUCT.LIMIT  1 
#define  PRINTOUTS  20.0 

(Function  declarations  23)  += 
void  initialize-parameters( ); 


based  upon  the  information  from  the  parsing  and 


40.  (Global  variables  4)  += 

real  sample-length  = 1.0;  /*  length  of  sample  in  meters  */ 

real  sulfate-failure.year\ 
real  chlonde-failure-year\ 

real  /c[NUM_CELLS];  /*  perm,  coefficient  due  to  cracking  */ 
real  sulfate-rate ; /*  dimensionless  rate  */ 

real  sulfate-raie-ms\  /*  m/sec  */ 

real  sulfate-depth  = 0.0;  /*  depth  of  sulfate  penetration  */ 

real  chloride-depth  = 0.0;  /*  depth  of  chloride  penetration  */ 

real  rebar-depth  = 1.0; 

boolean  termination  = FALSE;  /*  simulation  termination  */ 
int  termination-type  = TIME.LIMIT;  /*  limiting  condition  for  termination  */ 
MATERIAL  MaxDay  = {10000.0, TRUE};  /*  time  limit  (day)  */ 
real  MaxDepth  = 100.0;  /*  penetration  depth  limit  (m)  */ 

MATERIAL  crack.width  = {0.000100, TRUE}; 

real  crackspacing  = 100.0; 

real  crack-depth  = 0.0; 

real  crack-permeability; 

real  joint-width  = 0.0; 

real  joint-spacing  = 100.0; 

real  joint-lifetime  = 20000.0; 

boolean  joint-change-flag  = FALSE; 

boolean  joint-isspecified  = FALSE; 

real  joint-permeability  = 0.0; 

real  Pout\  /*  dimensionless  ext.  pressure  */ 

real  dP\  /*  dimensionless  rho  g h pressure  */ 

real  u[NUM_CELLS];  /*  Darcy  velocity  */ 

real  /i  = 0.001;  /*  pore  fluid  viscosity  */ 

real  g — 9.8;  /*  gravitational  constant(m/sec/sec)  */ 

real  p = 1000;  /*  density  of  pore  fluid  (kg/m3)  */ 

real  ks\  /*  bulk  permeability  */ 

real  Db\  /*  bulk  diffusivity  */ 
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real  <f>  = 0.16;  /*  porosity  */ 

real  F;  /*  formation  factor  */ 

real  a;  /*  degree  of  hydration  */ 

real  /;  /*  initial  solids  fraction  */ 

real  Vsample;  /*  volume  of  element  */ 

real  pch  = 33.1;  /*  density  of  CH  */ 

real  molesCH  = 13.0;  /*  moles  CH/ltr  soln.  */ 

real  i?o  = 0.00281;  /*  initial  1/F  */ 

real  r9;  /*  current  1/F  */ 

real  prini.day;  /*  next  day  to  print  interm,  results  */ 
real  prini.day  -interval]  /*  days  between  interm,  results  */ 

real  earliesi.failure.day ; /*  earliest  structural  failure  */ 


41.  initialize.parameters().  This  routine  initializes  the  scaling  parameters  and  physical  parame- 
ters. 

void  initialize.parameters( ) 

{ . 

int  i,  j,  km, 
int  n,  m; 
int  iter ; 

real  moles. cation , moles.anion , charge , x,  dx , z,  critical.Cl.concentration , T.cnt, 
concen.raiio , darcy .velocity , perm.time,  chloride.time,  perm.} actor , 

/*  modify  DT  to  account  for  high  permeability  */ 
perm-depth , ^4,  crack-factor ; /*  modify  DT  to  account  for  cracks  */ 

for  (fc  = FIRST_CELL;  k < NUH.SURFACES;  jb++) 

A X[k]  = 1.0; 

( Establish  material  parameters  of  concrete  42 ) 
aci.211  ( ); 

( Propagate  ion  concentration  in  concrete  45 ) 

(Print  parameters  to  stdout  47) 

(Initialize  ion  ‘eta’  parameter  48) 

L = sample-length / NUM_CELLS;  /*  the  universal  length  scale  */ 

MaxDepth  = MaxDepih  / L;  /*  dimensionless  */ 

(Calculate  crack  and  joint  adjustment  to  permeability  49) 
critical.Cl.concentration  = 0.000400*  Vsample  * (1.0  — </>)*  2.5/35.4; 
if  (anion  [C/].c[FIRST_CELL]  > critical.Cl.concentration)  { 

/*  calculate  T for  [Cl_crit]/[Cl(x=0)]  = erfc(X/2  sqrt(DT))  */ 
concen.raiio  = inv.erfc(critical.Cl.concentration  / anion  [C/].c[0]); 

T.crit  = SQR((rebar.depth / L)/(2.0  * concen.raiio))-, 
chlonde.failure.year  = 3.171  • 10~08  * T.crit  * SQR (L)/D0\ 

if  (head .value  > 0.0)  { /*  correct  diffusion  estimate  to  include  darcy  penetration  */ 

darcy.velociiy  = (k0/p)  * (p  * g * head  .value)/ sample-length; 

A = 1.0/ (4.0  * D0  * SQR  (concen.raiio)); 

perm.depih  = rebar.depth  + 1.0/ (2.0  * A * darcy.velociiy)  — sqrt(rebar.depth / (A  * 
darcy.velociiy)  + 1.0/ (4.0  * SQR(T  * darcy.velociiy))); 
perm.time  = perm.depih  / darcy.velociiy; 
chloride.failure.year  = 3.171  • 10-08  * perm.time; 

} 

} 

else  chlonde.failure.year  = 1.0-  10+1°; 
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print/ ("\"ChlorideulailureLJ(yr)\"\t\tu,/,8 . 011\n",  chloride-failure-year)-, 
perm-factor  = MAX(£o/3.0  • 10~18, 1.0)  * MAX(Aearf .value / 10.,  1.0); 
if  (crack.widih .value  = 0.0)  { 
crack-factor  = 1.0; 

} 

else  { 

if  ( sample-length  — crack-depth  > 0.05  * sample-length) 

crack-factor  = sample-length /{sample-length  — crack-depth ); 
else  crack-factor  = 10.0; 

} 

A T = 1. 0/(64. 0*  perm.factor  * crack-factor  * amon[Cl].r]  * i?o)i 

sulfate-rate-ms  = Youngs  .value  * SQB.(P. value)  * anion  [S04].c[FIRST_CELL]  * CE. value  * 

D0 /(roughness  .value  * 7. value  * (1  — v. value))',  /*  m/sec  */ 
sulfate-rate  = sulfate-rate.ms  * ( L/ Da)\  /*  AX/ AT  */ 

if  ( sulfate-rate  > 0.0)  sulfate-failure-year  = 3.171  • 10-08  * ( MaxDepth  * L)/ sulfate-rate-ms ; 
else  sulfate.failure-year  = 1.0-  10+1°; 

print/ Sul f ateuf  ailureu(yr  )\"\t\tu'/,8 . Olf \n",  sulfate.failure-year)-, 

/*  calculate  number  of  days  until  termination  */ 
earliest-failure-day  = 365.  * MIH  (sulfate-failure-year , chloride-failure-year)-, 

/*  determine  intervals  between  print  outs  */ 
print-day-interval  = Yll}l(earliest-failure-day , MaxDay  .value  )/PRINT_OUTS; 
print-day  = print-day-interval ; 

Pout  = ( k0/(p  * D0))  * p * g * head  .value-, 

dP  = (k0/(p  * D0))  * p * g * (sample-length /YIUH-CELLS)-,  /*  change  per  cell  */ 
for  (it  = 0;  k < NUM.SURFACES;  it++)  { 

m = do; 

<p'n[k]  = <j>-, 

<t>,n~l[k\  = <J)\ 

= Pout  * (real)  (NUM.SURFACES  - k - 1)/(NUM_SURFACES  - 1); 

} 

(Initial  estimate  of  pH  50 ) 

} 


42.  The  material  parameters  of  the  concrete  must  be  determined  based  upon  the  information  given 
by  the  user.  Missing  pieces  of  information  must  be  calculated  from  established  relationships.  The  . 
is-default  portion  of  the  material  data  determines  whether  the  user  specified  the  value. 

( Establish  material  parameters  of  concrete  42 ) = 

if  (Doo.is-de fault  / TRUE)  { /*  calculate  porosity  */ 

if  (wc .is-default  = TRUE)  wc. value  = (logl0(\0000.  * Doo.value)  + 9.84) / 6.0; 
if  (iboo .is-default  = TRUE) 

/*  kinfty. value  = 1.0E-18*(0.8904+.002525*exp(15.07*wc. value));  */ 
koo. value  = pow(10.0,  (5.0  * wc. value  — 21.0)); 

} 

if  (koo-is-default  / TRUE)  { 
if  (Doo. is-default  = TRUE)  { 

Doo.value  = 1.0  • IO-04  * pou>(10.0,  6.0  * wc. value  — 9.84); 

if  (wc.is-default  = TRUE)  wc. value  = (logl0(  10000.  * Doo.value)  -f  9.84)/ 6.0; 

} 

} 

if  (wc.is-default  / TRUE)  { /*  calculate  diffusivity  and  porosity  */ 
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if  {Doo.is.de fault  — TRUE)  .value  = 1.0  • 10  04  * pou;(10.0,  6.0  * wc. value  — 9.84); 
if  (k^-is-defaull  = TRUE) 

/*  kinfty. value  = 1.0E-18*(0.8904+.002525*exp(15.07*wc. value));  */ 
koo- value  = pou;(10.0,  (5.0  * wc. value  — 21.0)); 

} 

i?o  = Doq. value  / amon[Cl].Dj\ 

if  (d0  < DLIMIT)  d0  = DLIMIT  + 0.0001; 

if  (d0  < DLIMIT  + 0.07  * SQR(.18))  <t>  = sqrl((d0  - DLIMIT)/. 07); 
else  <p  = .17326  + sqrt(. 03002  - (.05832  + DLIMIT  - i?0)/1.87); 
a = (3.2  * wc. value)/ 1.36  — <f>  * (1  + 3.2  * wc. value)/ 1.36; 

/ = 1 .0/ ( 1 + 3.2  * wc. value)] 

Vsample  = 1000.0/</>; 

molesCH  = a * (0.61)  * f * Vsample / pen', 

F = 1.0/do; 

D0  — Doq.  value; 
k o = koo. value; 

/dre[0]  = 0.001  * Vsample  * <p; 
for  (k  = 1;  k < NUM.SURFACES;  jfc++)  { 
soLarray  [Ca]  [OH],  s[£]  = molesCH; 
litre[k\  = 0.001  * Vsample  * <p; 

} 

for  (k  = 0;  k < NUM_CELLS;  k++)  { 

£[*]  = ^o; 

<Pn[k)  = <P\ 

} 

This  code  is  used  in  section  41. 


43.  The  routine  aci.211  estimates  the  strength,  in  psi,  of  the  concrete  in  the  individual  com- 
putational elements.  The  estimate  is  based  upon  Table  5.3.4(a)  of  ACI  211.1-81  for  air-entrained 
concrete.  The  strength  can  be  adjusted  to  account  for  changes  in  porosity  using: 


pst'ft]  = psi0 


ln(<T[*1) 

In  (phi0) 


{ Function  declarations  23 ) += 
void  aci.211  (void); 


44.  void  aci-211{) 

{ m 

int  i; 
real  xp; 

xp  = pow  (exp  (wc.  value  — 3.575),  —2.6817); 

for  (t  = 0;  i < NUM.SURFACES;  *++)  sirength[i)  = t p; 

) 

45.  The  function  yyparse  sets  the  concentration  in  the  concrete  by  putting  the  user-specified  value 
into  the  . c[l]  cell.  The  internal  concentration  is  established  by  copying  this  value  into  the  remaining 
cells.  However,  this  procedure  is  not  performed  for  OH  and  H. 
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(Propagate  ion  concentration  in  concrete  45)  = 

for  (*  = 1;  i < num.cations]  i++)  { 

cation  [i].mo/es[0]  = cation  [*].c[0]  * /itre[0]; 
cation  [i].mo/es[l]  = cation\i\.c[  1]  * litre[\]\ 
for  (k  = 2;  k < NUM.SURFACES;  k++)  { 
cation[i].c[k\  = cation  [i].c[l]; 
cation[i\.moles[k ] = cation [z].c[fc]  * litre[k]\ 

} 

} 

for  (j  = 1;  j < num.anions ; j++)  { 

anion[j].moles[  0]  = anion  [j].c[0]  * /itre[0]; 
anion[j].moles[  1]  = anion  [;].c[l]  * /*"<re  [1]; 
for  (k  = 2;  k < NUH.SURFACES;  k++)  { 
anion  [j].c[fc]  = anion  [/].c[l]; 
anion[j].moles[k ] = anion[j].c[k]  * litre[k\] 

} 

} 

This  code  is  used  in  section  41. 


46.  Report  material  parameters  to  stdoui . 

(Preprocessor  definitions  11)  += 

#define  GET_STAT(a)  (( a.is.default  = TRUE)  ? "DEFAULT"  : "USER") 


47.  (Print  parameters  to  stdoui  47)  = 

prini/(  "\"THICKNESSu\"\tu'/.8 . 51f  \t\"  (m)uuuuuu\"\t\"'/,s\"\n",  sampleJcngth , 

GET_STAT(  thickness))] 

print/(»\"DIFFUUUUUu\"\tu7.8.11e\t\"(in-2/sec)\"\t\",/.s\"\n",DCJ,GET_STAT(Doo)); 
pnntf  {"\" PERMuuuuuu\"\tu,/.8.11e\t\"(m/sec)uu\"\t\",/,s\"\n",9.8  • 10+o6  * k0, 
GET_STAT(&oo)); 

print/ ("\"WCUUUuuUUu\"\tu,/*8 . 51f  \t\"uuuuuuuuu\"\t\",/,s\"\n",  wc. value , GET_STAT(u>c)); 
print/ ("\"HEADLJUUUUL|\"\tu,/,8 . 51f  \t\"  (m)uuuuuu\"\t\"y,s\"\n",  head  .value , GET_STAT(/iead)); 
print/  ("Snli  ateuAttackuParameters : \n"); 

print/ ("V'  Y0UNGSuuuu\"\tu,/.8.11e\t\"(N/m‘2)uu\"\t\",/.s\"\n")  Youngs  .value , 

GET_STAT(  Youngs))] 

prin//("\"BETAUULJLjLjU\"\tu'/,8.  lle\t\"uuuuuuuuu\"\t\"*/,s\"\n",  /?.va/ue , GET_STAT(/?)); 
print/("\"CEuuuuuuuu\"\tu'/.8 . 51f  \t\"  (Mol/m*3)\"\t\"*/,s\"\n",  CE.  value , GET_STAT(CE)); 
print/ ("\"R0UGHNESSu\"\tu'/,8 . 51f  \t\"uuuuuuuuu\"\t\"'/,s\"\n",  roughness  .value , 

GET _ST kT(roughness )); 

pnnt/("\"GAMMAuuuuu\"\tu'/,8 . 511 \t\"  ( J/m~2)uu\"\tV7.s\"\n",  y. value , GET_STAT(t)); 
prin</( "\ "P0ISS0Nuuu\"\tu'/,8. 51f\t\"uuuuuuuuu\"\t\"'/,s\"\n",  v. value,  GET_STAT(i/)); 
pnnt/l"\ n"); 

if  ( MaxDepth  < sample-length ) 

pnnt/("\"DEPTHULILJUU\"\tu,/.8 . 51f \t\"  (m)  uuuuuu\"\t\",/.s\"\n" , MaxDepth , "USER"); 
if  ( rebar.depth  < sample-length) 

print/ ("\"REBARUUUUu\"\tu'/.8 . 51f \t\"  (m)UUUUuuV7tV7.s\"\n" , rebar.depth , "USER"); 
print/ ("\"TIMEuuuuljLj\"\tu,/,8 . Olf  \t\"  (day  )uuuu\"\t\"'/,s\"\n",  MaxDay  .value , 

GET  _STkT(  Max  Day )); 

if  (crack-width  .is.de/ault  = TRUE  A neutraLaxis.depth  .is.de/ault  = TRUE)  { 
crack-depth  = 0.25  * sample-length] 
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crackspacing  = 2.0; 

} 

if  ( neutraLaxis.depih.is.de/ault  = FALSE) 

crack-depth  = sample-length  — neuiraLaxis-depih  .value ; 
if  (crack-width. value  >0.0)  prz7i//(,,\"CRACKu=u,/.8.511uATu,/.8.511uDEPTHu,/.8.51f\M\nH, 
crack-width  .value , crack-spacing , crack-depth ); 
if  ( joint-permeability  > 0)  { 

pnntf("\"  J0INTuPERMu=u'/,8 . lle\"\n" , joint-permeability)-, 

pnntf("\"  J0INT='/,8 . 51f  uATu*/,8 . 51euUNTILu’/.8 . 511  \"\nM , joint-width,  joint-spacing , 
joint-lifetime /365); 

} 

This  code  is  used  in  section  41. 


48.  Once  the  parameter  D0  has  been  established,  the  value  of  ij  for  each  ion  must  be  initialized. 
(Initialize  ion  ‘eta’  parameter  48)  = 
for  (i  = 0;  i < num.cations;  z++)  cation[i].rj  = cation[i\.Dj / D0\ 
for  (j  = 0;  j < num.anions\  j++)  anion  [j]. tj  = anion  [j].Dj / D0\ 

This  code  is  used  in  section  41. 


49.  The  presence  of  cracks  is  reflected  in  k which  a multiplicative  adjustment  for  permeability. 
In  the  absence  of  cracks  «[i]=l.  The  quantities  crack.width  and  crackspacing  are  specified  by  the 
user. 

The  equation  for  the  permeability  of  cracks  of  width  w spaced  a distance  a apart  is 


(Calculate  crack  and  joint  adjustment  to  permeability  49)  = 
crack-permeability  = (SQR(crack-width  .value) / 12); 
x = 0;  /*  x starts  from  inside  surface  */ 

dx  = sample-length  / NUM_CELLS; 
for  ( k = NUH.CELLS  - 1;  k > 0;  k — ) { 
x +=  dx\ 

if  (x  < crack-depth)  { 

k[&]  = k0  * (1.  — (joint-width / joint-spacing)  — (crack.width .value / crackspacing))] 

/*  weighted  length  */ 

k[&]  +=  joint-permeability  * (joint-width /joint-spacing)-, 

-f=  crack-permeability  * (crack.width .value / crackspacing)] 

} 

else  { 

if  (x  — dx  < crack-depth)  { /*  crack  ends  in  this  cell*/ 

n[k]  = k0  * (1.  — (joint-width /joint-spacing)  — (crack.width. value / crackspacing))] 

/*  weighted  length  */ 

/c[Ar]  +=  joint-permeability  * (joint-width / joint-spacing)] 
k[&]  = (x  — crack-depth)  / n[k]\ 

k[&]  +=  (dx  — (x  — crack-depth)) / (crack-permeability  * (crack.width. value / crackspacing)) 
k[&]  /=  dx\ 
k[&]  = 1.0/#c[Ar] ; 

} 

else  { 


' 


' 
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K,[k]  = k0  * (1.  — (joint-width  / joint-spacing)); 

k[£]  +—  joint-permeability  * (joint-width /joint-spacing); 

} 

} 

K[k]  = K[k]/k0; 

} 

This  code  is  used  in  section  41. 


50.  Need  to  make  an  initial  estimate  of  the  equilibrium  concentrations  and  the  pH  of  each  cell. 
(Initial  estimate  of  pH  50 ) = 

/*  make  a rough  guess  for  concentrations  in  cases  of  solids  present  */ 
for  (k  = 0;  k < NUM.SURFACES;  k++)  { 
for  (i  = 0;  i < num-cations;  i++)  { 
for  (j  = 0 ; j < num-anions;  j++)  { 
if  (so/_arraj/[i][j]..s[&]  > 0.0)  { 
m = sol-array[i][j].m; 
n = soLarray  [i)[j).n; 
if  (m  = 0 V n = 0)  { 

if  (m  = 0)  pnntf("m=Ou\rx"); 
if  (n  = 0)  pnntf(" n=Ou\n"); 

pnnt/("\n\tProblemuwithu*/,s*/,s\n" , cation[i\.name , anion  [j].name); 
exit(  0); 

} 

moles-cahon  = pow  (soLarray  [i]\j]. ksp , 

1./ (m  4-  n))  * pow  (anion  [j],  valence  / caiion[i\.  valence , n/(n  + m)); 
moles.anion  = moles-cation  * (cation[i\. valence / anion[j). valence); 
cation[i\.c[k\  = moles.cation; 
anion[j].c[k]  = moles.anion; 
cation[i].moles[k]  = cation[i].c[k]  * htre[k\; 
anion[j].moles[k]  = anion  [i].c[fc]  * litre[k]; 

} 

} 

} 

} /*  calculated  [OH]  and  [H]  concentrations  */ 

for  (k  = 0;  k < NUM_SURFACES;  k++)  { 
charge  = 0.0; 

for  (j  = 1;  j < num.amons;  ji++)  charge  +=  anion[j]. valence  * anion\j\.c[k\; 
for  (i  =1;  i < num.cations;  i++)  charge  —=  cation[i\. valence  * cation[i\.c[k\; 
anion  [OH] . c[fc]  = 0.5  * ( — charge  + sqri(SQK(charge)  + 4.0  * soLarray  [//][0H].£sp)); 
cation[H].c[k]  = soLarray  [//][OH].A:sp/am0n  [OH]  .c[Ar] ; 
anion  [OH]. moles[k]  = anion  [OH] .c[Ar]  * htre[k\; 
cation[H].moles[k]  = caiion[H].c[k]  * litre[k\; 

} 

This  code  is  used  in  section  41. 


51.  Inverse  complementary  error  function.  Find  x given  C such  that  C = erfc(x). 

(Function  declarations  23)  += 
real  inv.erfc( real); 
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52.  Use  lookup  table  from  data  in  Table  7.1  from  Abramowitz  and  Stegun,  Handbook  of  Math- 
ematical Functions.  Get  bounds  on  y data  and  linear  interpolate  x data.  The  variable  frac  is  the 
fraction  of  the  gap  between  neighboring  y values.  Note:  additions  to  the  table  require  changes  to 
the  main  routine  to  reflect  the  new  dimension  of  the  arrays, 
real  tnv.erfc{ real  C) 

{ 

real  x,  frac , x.iable[ ) = {0.00,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90, 1.00, 1.20, 1.40, 
1.60, 1.80,  2.00},  y.lable[]  - {1.00,  .888,  .777,  .671,  .572,  .480,  .396,  .322,  .258,  .203,  .157, 

.090,  .048,  .024,  .011,  .005}; 

int  t; 

i = 16; 

while  (t  > 0 A y.iable[ — i]  < C)  ; /*  find  bounds  */ 

if  (x  < 15)  { 

frac  = ( y.iable[i ] — C)/(y.table[i\  — y.iable[i  + 1]); 
x = x.iable[i]  + frac  * ( x.table[i  + 1]  — x.table[i})\ 

} 

else  x = x_fa6/e[15];  /*  if  C > 2.0  simply  use  last  x value  */ 

return  x; 

} 
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53.  ADVECTION- DIFFUSION. 

Ion  Transport 

At  the  core  of  4SIGHT  is  the  advection-diffusion  equation  to  account  for  both  diffusion  of  ions 
and  Darcy  flow  of  the  pore  solution  due  to  hydrostatic  head.  The  flux  of  ions  due  to  both  gradients 
in  the  ion  concentration  and  to  a volume  average  flow  of  pore  solution  is 

j = — DVc  -f-  cu 

where  j is  the  ion  flux,  D is  the  diffusivity,  c is  the  ion  concentration,  and  u is  the  volume-averaged 
velocity  of  the  pore  solution.  The  time  dependent  change  in  concentration  is  the  negative  divergence 
of  the  flux: 

dc 

— = V ■ DVc  — u • Vc  — cV  • u 
dt 

Given  a hydrostatic  pressure  head  on  a vertical  column  of  porous  media,  the  pore  volume-averaged 
flow  V£>  is 

vD  = (Vp  - pg) 

For  the  hydrostatic  heads  considered  here,  the  body  force  term,  pg  is  non-negligible.  This  equation 
can  be  cast  into  the  more  familiar  form,  assuming  a constant  density  pore  fluid,  using  a modified 
pressure  potential: 

^ = p - pgz 

This  gives  the  more  familiar  Darcy  equation 


\D  = 

F 

The  volume-averaged  velocity  u can  be  related  to  the  Darcy  flow  velocity: 

V D = 


where  <p  is  the  porosity. 

Finally,  the  above  equations  can  be  combined  to  give 

c)c  h Ic 

— = V ■ DVc  + — Vip  • Vc  + cV  ■ — Vip 
dt  <pp  <pp 

This  equation  gives  the  spatial  and  temporal  behavior  of  the  concentrations.  To  complete  the 
calculations  a means  is  needed  to  update  the  hydrostatic  pressure  potential,  ip. 

Continuity  Equation 

The  temporal  behavior  of  ip  is  calculated  using  the  continuity  equation: 


dp 

dt 


= —V  ■ p\ 


where  v is  the  intrinsic  velocity  of  the  pore  solution.  After  averaging  over  the  microstructure,  the 
continuity  equation  becomes 

d<P  _ 

a=-Vvo 

This  can  be  related  back  to  pressure  using  the  Darcy  equation  once  again: 


-Vip 
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54.  Ion  transport. 

(Function  declarations  23)  += 
void  ion-diffusion  (void); 


55.  ion-diffusion. 
void  ion-diffusion  () 

{ 

int  *,  j,  k\ 

real  cation-flux [10][mm_S\JRFkCES],  anion.flux [10][NUM_SURFACES],  upper-limit  = 5.0  • 10+o\ 
lower-limit  = 5.0  • 10+o3,  max-flux-ratio\ 
real  tmp.cation  [10][NUM_SURFACES],  tmp.amon  [10][NUM_SURFACES]; 

AX[FIRST_CELL]  = (real)  (FIRST.CELL  + 1)  - sulfate-depth ; 

max -flux  .ratio  =■  0.0; 

for  (i  = 0;  i < num-cations\  i++)  { 

for  ( k = FIRST.CELL  + 1;  k < NUM.SURFACES  - 1;  k++)  { 

cation-flux  [z][fc]  = cation[i].r]  * (£[&]  * ( cation  [i].c[fc  + 1]  — cation  [i].c[£])/AX[fc]  — f[k  — 1]  * 
( cation [i].c[^]  — cation[i].c[k  — 1])/AX[/:  — 1])/ (0.5  * (AX[£]  + AX[fc  — 1])); 
cation-flux  [i\[k]  —=  (1  ./<f>,n[k])  * * (( cation[i].c[k  + 1]  — cation[i].c[k  — l])/(A>Y[fc]  + 

AX[t-l])); 
if  ( cation [i].c[Ar]  > 0.0) 

max-flux-ratio  — MkX(max-flux-ratio , AT  * cation-flux  [i][k\/ cation[i\.c[k])\ 

} 

} 

for  ( j = 0;  j < num.amons ; j++)  { 

for  ( k = FIRST.CELL  + 1;  k < NUM.SURFACES  - 1;  k++)  { 

anion.flux  [j][k]  = anion  [ j\.i)  * (£[&]  * ( anion[j\.c[k  + 1]  — anion[j].c[k])/ AX[k]  — £[k  — 1]  * 
(anion [j].c[fc]  — anion[j].c[k  — l])/A.Y[fc  — 1])/ (0.5  * (AX[fc]  + AX[1:  — 1])); 
anion-flux  [_/] [A:]  — = (1  ,/(f)'n[k])  * t>[&]  * ((anion  [ji].c[fc  + 1]  — anion  [?‘].c[fc  — 1])/(A.Y[&]  + 
AX[i-l]))i 
if  (anion  [i].c[fc]  > 0.0) 

max-flux-raiio  = Vlhl(max-flux-ratio , AT  * anion-flux  [j][k\/ anion[j].c[k])\ 

} 

} 

if  (max-flux-raiio  > upper-limit)  AT  *=  upper-limit / max-flux.ratio; 
if  ( max-flux-raiio  < lower-limit ) AT  *=  upper-limit / max-flux-raiio \ 
for  (k  = FIRST.CELL  + 1;  k < NUM_SURFACES  - 1;  k++)  { 

for  (i  = 0;  i < num-cations ; i++)  cation[i].c[k\  +=  AT  * cation-flux [i][k]\ 
for  (j  = 0;  j < num.anions\  j++)  anion[j].c[k]  -f=  AT  * anion-flux [j][k]\ 

} 

k = NUM.SURFACES  - 1; 

for  (i  = 0;  i < num-caiions\  i++)  cation[i].c[k)  = cation[i].c[k  — 1]; 
for  (j  = 0;  j < num-anions\  j++)  anion[j].c[k]  = anion[j].c[k  — 1]; 
for  (i  = 0;  i < num.caiions\  i++) 

for  (k  = FIRST.CELL;  k < NUM.SURFACES;  k++) 
cation[i].moles[k]  = cation[i\.c[k\  * htre[k]  * AY[ij; 
for  (j  ~ 0 ; j < num-anions ; j++) 

for  (k  = FIRST.CELL  + 1;  k < NUM.SURFACES;  k++) 
anion[j].moles[k]  = anion[j].c[k\  * litre[k]  * AY[fc]; 
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56.  Chemical  equilibrium. 

Given  the  number  density  of  ions  in  a computational  element,  determine  if  any  of  the  ions  should 
go  in  to/out  of  solution.  Determine  the  pH  and  adjust  H and  OH,  accordingly. 

(Function  declarations  23)  += 

void  chemical-equilibrium  (boolean); 


57.  Minimizing  function.  This  is  the  function  to  minimize  for  the  chemical-equilibrium  routine. 

The  objective  of  chemical-equilibrium  is  to  determine  how  many  moles  of  salt  should  be  leached/precipitated. 
Therefore,  minfunc  must  adjust  for  the  pore  volume. 

The  equation  minfunc  is  minimizing  the  square  of 

[c+™  ]"V  + ^]n  — *.,  = <> 

where  x is  moles  salt,  m and  n are  stoichiometric  ratios,  and  V is  the  pore  solution  volume, 
real  min/u7ic(real,  int,  int,  int); 
real  minfunc( real  x,int  i,int  j,int  k) 

{ 

real  C,  A; 
int  m,  n; 

m = soLarray[i][j].m\ 
n = soLarray[i][j].n\ 

C = cation[i].c[k]-, 

A = anion  [j].c[k]-, 

return  SQR(potr(C  + m*  x/(litre[k\  * AA'[fc]),m)  * pow(A  + n * x/(litre[k\  * AA'ffc]), 
n)  — sol-array  [i][j].ksp)\ 

} 


58.  chemicaLequilibnum( ).  Cycle  through  sol.array  and  determine  if  any  ions  should  go  in  to/out 
of  solution  based  upon  the  concentration,  the  solubility  constant,  and  the  presence  of  solid  salt. 

#defineTOL  1.0  • 10~O6 

void  chemical-equilibrium  (boolean  Change-Porosity) 

{ > 

int  i,  j,  k,  n,  m,  iterations ; 
real  w.max , <5; 

real  xa,  xb,  xc,  fa,  fb,  fc,  moles,  tmp\ 
real  dsolid ; 

real  A litre;  /*  change  in  pore  volume  */ 

real  old-cOH\ 
real  charge , temp-, 

boolean  NEG_FLAG; 

for  (k  = FIRST_CELL  + 1;  k < NUH.SURFACES;  k++)  { 
iterations  = 0; 
do  { 

W-max  = 0.0; 

for  (j  = 0;  j < num-amons ; j++)  { 
for  (i  = 0;  i < num-cations-,  i++)  { 
if  (~i(i  = 0 A j = 0))  { 
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n = so/_arraj/[z][j].n; 
m = soLarray[i][j].m ; 

temp  = pow(cation  [i].c[fc],  (double)  m)  * pow(anion  [jr] .c[Ar] , (double)  n); 
if  ((so/_arrap[z][/].s[&]  > 0.0)  V ( temp  > soLarray[i][j].ksp))  { 
xa  = — MIN(c<rfton[t].c[fc],  anion[j].c[fc])/10.0; 
xb  = —xa ; 
xc  = 0.0; 

mnbrak(Lxa , k.xb , kxc , &:/a , k.fb , hfc , mnifunc,  i,j , fc); 

/mp  = brent  (xa,xb,xc,minfunc,TQL,&.moles,i,j,k); 
if  (fabs(moles)  > 0.05) 

mo/es  *=  0.10;  /*  no  drastic  changes  */ 

if  (so/_array[z][/].s[fc]  — moles  < 0.0) 

moles  = so/_array  [*][;]  .s[fc];  /*  insufficient  salt  */ 

if  ( cation[i\.moles[k ] + m * moles  < 0.0) 
moles  = — cation[i].moles[k]/m  + 0.000001; 
if  ( anion [j].mo/es[fc]  + n*  moles  < 0.0) 
moles  = — anion  [j].moles[k]/n  + 0.000001; 

A litre  = 0.001  * moles  * sol.array[i][j].molar.density  / AX[k]; 
if  (litre[k]  + A litre  < 0.0)  { 

A litre  = — litre[k];  /*  insufficient  pore  space  */ 

moles  = Alitre  * AX[&]  * 1000.0/ sol.array[i\[j].molar.density; 

} 

cation[i\.moles[k } -f=  m*  moles ; 
anion\j].moles[k ] +=  n*  moles; 
if  (Change. Porosity  = TRUE)  { 
htre[k]  +—  Alitre ; 
sol.array  [z]  [/]  ,s[k]  — = moles ; 

} 

6 = fabs(moles); 
w.max  = MAX(<5,  w.max); 

) 

} 

} 

} /*  readjust  the  concentrations  */ 

4>'n[k]  = 1000.0  * hire  [A:]  * AX [k\/  Vsample; 

litre[k]  = MAX(/zire[fc],  0.0001);  /*  avoid  zero  volume  problems  */ 

for  (i  = 0;  i < num.caiions ; i++) 

cation[i].c[k]  = caiion[i\.moles[k\/(litre[k]  * AAr[£]); 
for  ( j = 0;  j < num.amons ; j++) 

amon\j].c[k\  = amon[j].moles[k]/(litre[k\  * AX[£j); 

/*  Determine  [OH]  and  [/^concentration  */  /*  store  old  value  */ 

old.cOH  = anion  [OH]. c[Ar] ; 
charge  = 0.0; 

for  (j  = 1;  j < num.amons;  j++)  charge  +=  anion  [j], valence  * anion[j].c[k\; 
for  ( i = 1;  i < num.caiions;  i++)  charge  —=  cation[i\. valence  * cation[i].c[k\; 
anzon[0H].c[/:]  = 0.5  * ( — charge  + sqri(SQR(charge)  + 4.0  * sol.array[H][OR].ksp)); 
anion[OH].mo/es[fc]  = anion  [OH].  c[fc]  * litre[k]  * AX[&]; 
if  (fabs((anion  [OH] .c[Ar]  — old.cOH)/ old.cOH)  > 1.00)  { 
w.max  = 1.0; 

anion  [OH] . c[fc]  = 0.5  * (anion  [OH], c[fc]  + old.cOH); 
anion  [OH]. moles[k]  = anion  [OHj.cffc]  * litre[k]  * AX[fc]; 
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} 

cahon[H].c[k]  = soLarray[H][QH].ksp  / anion  [OH].c[fc]; 
cation[H].moles[k ] = cahon[H].c[k]  * htre[k]  * AX  [A;]; 

} while  (( w.max  > 0.0001  V anion  [OH].c[A:]  < 0.0)  A iterations  ++  < 50); 

} /*  extrapolate  <?!>,n[FIRST_CELL]  value  */ 

<j)'n [FIRST_CELL]  = ^'"[FIRST.CELL  + 1]  - AX[FIRST_CELL]  * (<£'n[FIRST_CELL  + 2]  - 
</>m[FIRST_CELL  + 1]); 

<£'n[FIRST_CELL]  = MAX(i'n[FIRST_CELL],  0.0);  /*  extrapolate  <i>'n [NUM.CELL]  value  */ 

<j>'n [NUM_CELLS]  = 2.0  * <?i,n[NUM_CELLS  - 1]  - <?i/n[NUM_CELLS  - 2]; 

/*  extrapolate  /dre[FIRST_CELL]  value  */ 
litre [FIRST_CELL]  = 0.001  * Vsample  * <£'n[FIRST_CELL]; 

} 


59.  Minimization.  The  routine  mnbrak  brackets  a minimum. 

(Function  declarations  23)  += 

void  mnbrak( real  *,real  *,real  *,real  *,real  *,real  *, real(*/unc)(real, int, int, int), int, 
int,  int); 

^include  "mnbrak.  c" 


60.  Parabolic  approximation.  The  routine  brent  uses  parabolic  approximation  to  determine  the 
minimum  given  three  points  bracketing  the  minimum. 

(Function  declarations  23)  += 

real  brent  (real,  real,  real,  real(*/unc)(real,  int,  int,  int),  real,  real  *,  int,  int,  int); 

^include  "brent. c" 
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61.  Advance  sulfate  front. 

(Function  declarations  23)  += 
void  advance.sulfaie-front()\ 


62.  void  advancesulfate.front( ) 

{ . 

int  i,  j,  k,  new.first.cell ; 
sulfate-depth  +=  sulfate.rate  * AT; 

if  ( sulfate-depth  > (real)  (FIRST.CELL  + 1)  V AX[FIRST_CELL]  < 0.10)  { 

new-first-cel l = (int)  ( sulfate-depth  +0.4);  /*  penetration  beyond  next  cell  */ 

/*  copy  info  from  old  first  cell  into  new  first  cell  */ 
for  ( k = FIRST_CELL  + 1;  k < new.first.cel  l ; k++)  { 

for  (i  = 0;  i < num.cations ; i++)  cation[i\.c[k]  — cafion[f].c[FIRST_CELL]; 
for  ( j = 0;  j < num. unions',  j++)  anion[j].c[k]  = anion  [j].c[FIRST_CELL]; 
for  (i  = 0;  i < nuni-cations ; i++) 

for  (j  = 0;  j < num.amons',  j++)  so/_arTa2/[i][ji].s[fc]  = soLarray[i\ [j].s[FIRST_CELL]; 

} 

^[new.first.cell]  = ^[FIRST.CELL]; 

FIRST.CELL  = new.first.cell ; 

AA[FIRST_CELL]  = (real)  (FIRST.CELL+  1)  - sulfate.depth ; 

} 

} 
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63.  Chloride  penetration. 

The  ratio  of  the  mass  of  chloride  ions  to  the  mass  of  concrete  equal  to  0.0004  is  the  depth  of  chloride 
penetration.  The  ration  equals  10_3[C7~]Ac/<?V((l  — </>)2.5).  Aci{gm/mole)  is  the  gram  atomic 
mass  of  chlorine,  <t>  is  porosity,  and  2.5  is  the  density  of  concrete. 

( Global  variables  4 ) += 

real  CLralio  [NUM.SURFACES]; 


64.  (Assess  chloride  penetration  64)  = 

for  (k  = FIRST_CELL;  k < NUM.SURFACES;  k++) 

CLratio[k]  = 35.5  * anion[Cl].moles[k]/(2.b  * (1.0  — </>/n[£])  * Vsample); 
k - NUM.SURFACES  - 1; 
while  (k  > 0 A Cl.ratio[k]  < 0.0004)  k — ; 
if  ( k < 0)  chloride. depth  = 0.0;  /*  k at  lower  limit  */ 

else  { 

if  ( k = FIRST_CELL)  { /*  account  for  distance  equals  AX  */ 

chlonde.depih  = AX[FIRST_CELL]  * (CLratio[k]  — 0.0004)/(C7_raho[fc]  — CLratio[k  + 1]); 
chlonde.depih  +=  (real)  FIRST_CELL; 

} 

else  chlonde.depih  = (real)  k + ( Cl.ratio[k\  — 0.0004)/(C/_raho[A']  — Cl.raiio[k  + 1] ) ; 

} 

This  code  is  used  in  section  10. 
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65.  Update  Pressures. 

The  pressures  are  updated  using  the  continuity  equation: 

dp 

- = -V.„v 


where  v is  the  intrinisic  pore  fluid  velocity  and  p is  the  pore  fluid  density.  An  equation  for  the  bulk 
material  can  be  obtained  from  a volume  average  over  a represetative  volume  V : 


_1_ 

V 


L 


dt  V Jv 


V • p\  d3x 


For  these  equations,  assume  p is  a constant  since  water  is  virtually  incompressible.  Rearranging  and 
simplifying  the  above  equation  gives: 


d_  l 
dt  V 


-V 


U 


px  d3x 


Since  p is  zero  outside  the  pore  volume,  if  Vp  represents  only  the  pore  volume  then  the  above  equation 
simplifies  to 

d<i>  <t>  [ , 3 

dt  Vp  Jv 


which  finally  gives 


d<b 

— = -V  • 4>u  = -V  • xD 
dt 


where  xp  is  the  Darcy  velocity.  Substituting  for  the  Darcy  velocity  gives 

d(f)  — k 
— zz  V • -Vp 
dt  p 


( Function  declarations  23 ) += 
void  updaie.pressures  (void); 


66.  update.pressures . 
void  update.pressures{) 

{ _ 

int  k,  iteration ; 

real  perm  [NUM_CELLS],  imp  [NUM.SURFACES],  tol  = 0.0005,  e; 
for  (k  = FIRST.CELL;  k < NUM.CELLS;  k++)  perm[k]  = «[Jfe]  * CUB(£[Jfe]/tf0); 
[FIRST.CELL]  = Pout  + (NUM.CELLS  - 1 + A A [FIRST.CELL])  * dP\ 
[NUM.SURFACES  - 1]  = 0; 
iteration  = 0; 
do  { 

for  {k  = FIRST.CELL  + 1;  k < NUM.SURFACES  - 1;  Jb++)  { 

imp  [A:]  = perm  [fc]  * ^[k  -f  1]/AA[&]  4-  perm  [fc  — 1]  * \t[A  — 1]/AX[&  — 1]; 

/*  ****  tmp[k]  -=  0.5*(DX[k-l]+DX[k])*(PHIn[k]-PHIn[k-l])/DT;  ***  */ 
tmp[k ] *=  1.0/(perm[A]/AX[A]  4-  perm[k  — l]/AX[fc  — 1]); 

} 

e = 0.0; 

for  ( k = FIRST.CELL  + 1;  k < NUM.SURFACES  - 1;  k-H-)  { 
if  [A:]  > 0.0)  e = MAX(e,/a6s(^[/:]  - tmp  [A]) /$[/:]); 
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'!'[/:]  = tmp[k]\ 

} 

} while  (e  > iol  A iteration ++  < 1000); 
for  (k  = 0;  k < NUM_CELLS;  k++) 

f [fc]  = -K[k]  * CUB(^[jb]/i9o)  * (VI >[k  + 1]  - *[*])/ AX[*]; 


67.  Copy  the  new  porosities  into  the  n — 1 values. 

(Update  the  porosities  67)  = 

for  (Jfc  = FIRST.CELL;  k < NUM.SURFACES;  k++)  <£m_1[fc]  = 4'n[k] ; 
This  code  is  used  in  section  9. 
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68.  Adjust  physical  parameters. 

Due  to  dissolution/precipitation. 

(Adjust  physical  parameters  68)  = 

for  ( k = FIRST.CELL;  k < NUM_CELLS;  k++)  { 

</>"[£]  = 0.5  * (<^,n[fc]  + <t>'n[k  + 1]);  /*  interpolate  phi’  values  */ 

i?  = DLIMIT  + 0.07  * SQR(d>n[/;]); 
if  (<j>n[k]  > 0.180)  d +=  1.8  * SQR (<j>n[k)  - .180); 
if  (<£"[&]  > <f>)  /*  leaching  */ 

£[*]  = + 5.0  * (tf  - tf0); 

else  /*  precipitation  */ 

£[*]  = 

} 

£'[FIRST_CELL]  = £[FIRST_CELL]  - (AX[FIRST_CELL]/(1.0  + AA[FIRST_CELL]))  * (£[FIRST_CELL  + 
1]  - £[FIRST_CELL]); 

for  ( k = FIRST.CELL  +1;  k < NUM_SURFACES  - 1;  k++)  = 0.5  * (£[k]  + £[Jfe  - 1]); 

£'[NUM_SURFACES  - 1]  = 1.5  * £[NUM_CELLS  - 1]  - 0.5  * £[NUM_CELLS  - 2]; 

This  code  is  used  in  section  9. 


' 
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69.  Proper  time. 

The  clock  time  must  be  advanced 
(Global  variables  4)  += 

real  Time  = 0;  /*  dimensionless  cumulative  time  */ 

real  Day  = 0;  /*  cumulative  time  (day)  */ 


70.  (Advance  global  clock  7o)  = 
Time  +=  AT; 

Day  = Time  * SQR (L)/(D0  * 86400); 

This  code  is  used  in  section  9. 


71.  Check  time  dependent  variables. 

(Check  time  dependencies  71 ) = 

if  (Day  > joint-lifetime  A joint.is.specified  = TRUE)  { 
joint.change.flag  = TRUE; 

} 

This  code  is  used  in  section  10. 


72.  ( Assess  simulation  termination  72  ) = 

if  (Day  > MaxD ay  .value)  termination  = TRUE; 
if  (sulfate. depth  > MaxDepth ) termination  — TRUE; 
if  ((chloride-depth  * L)  > rebar.depth)  termination  = TRUE; 
This  code  is  used  in  section  10. 
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73.  OUTPUT. 

The  following  are  output  routines. 


74.  Print  the  state  of  the  system  for  the  first  4 concentraction  surfaces. 

(Print  initial  system  state  74)  = 
prmf/("\n\nInitialustateuoflJsystem:  \n\n"); 
pn  7if/("\"uI0Nu\"\t\"EXTERNAL\"\t\"  INTERN  AL\" \n" ); 
for  (i  = 0;  i < num.caiions\  i++)  { 

if  ( cation [i].c[FIRST_CELL]  ^ 0 V cation [i].c[FIRST_CELL  + 1]  ^ 0)  { 
pnntf  ("\"'/,4s cation[i].name)\ 

for  (k  = FIRST.CELL;  k < FIRST_CELL  + 2;  Jfc++)  print/ ("7.8 . 51f\t",  cation[i\.c[k])\ 
print/ ("\ n"); 

} 

} 

for  (j  = 0;  j < num.amons ; j++)  { 

if  (anion  [/].c[FIRST_CELL]  ^ 0 V anion  [j].c[FIRST_CELL  + 1]  ^ 0)  { 
prjnt/("\"'/,4s : anion[j].name)\ 

for  (k  = FIRST_CELL;  k < FIRST_CELL  -f  2;  k++)  print/ ("7.8 , anion[j].c[k]); 
print/ ("\ n"); 

} 

} 

pnnf/("\"uupH:\"\tM); 

for  (k  = FIRST_CELL;  k < FIRST_CELL  4-  2;  k++)  print/ ("7.8 .51i\t" , -loglO (cation[H].c[k}))\ 
pnnt/("\  n"); 
print/  ("\n"); 

This  code  is  used  in  section  7. 


75.  Print  header  for  permeability,  sulfate  penetration,  and  chloride  penetration  depths. 

( Print  header  75 ) = 

pr2nf/("\"Day\"\t\"L\"\t\"Kuu\"\t\"DUUUu\"\t\"S04\"\t\"Cl\"\t\"Flux\"\tuu\"pH\"\n"); 
print/ ("\"uuu\"\t\"m\"\t\"m/ s\"\t\"m-2/ s\"\t\"muu\"\t\"mu\"\t\"ml/dy/m2\"\n"); 
print/  ("\n"); 

This  code  is  used  in  section  2. 


76.  Print  intermediate  results  for  permeability,  sulfate  and  chloride  penetration  depths.  The  bulk 
values  of  permeability  and  diffusivity  must  be  calculated  here. 

(Print  intermediate  results  76)  = 
if  ( Day  = 0.0  V Day  > print-day ) { 
print-day  +=  pnnt-day-interval\ 

print/ ("7. ldu\t'/.6 . 31f  ",  (long)  Day  , (AX[FIRST_CELL]  + NUM_CELLS  - FIRST.CELL  - 1)  * L)\ 
kB  = AX[FIRST_CELL]  * CUB (i?o/£ [FIRST_CELL] ) / ac[FIRST_CELL] ; 

Db  = AX  [FIRST.CELL]  * (t?0/£[FIRST_CELL]); 
for  (k  = FIRST.CELL  + 1;  k < NUM.CELLS;  k++)  { 

kB  +=  cub(iV£[*])M*]; 
db  +=  0o/£[*]; 

} 

kB  = (AX[FIRST_CELL]  + NUM.CELLS  - FIRST.CELL  - 1)  * k0/kB\ 
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Db  = (AX[FIRST_CELL]  + NUM.CELLS  - FIRST.CELL  - 1)  * Dd/Db; 
print/ ("\t7,6 . lie",  9.8  ■ 10+O6  * kB); 
pnntf  ("\t7,6 . lie" , DB); 

print/ ("\t7,6 . 31f  \t'/,6 . 31f  ",  sulfate.depth  * L,  chloride-depth  * L); 
print/ ("\t'/,6 . 31f  " , 1000.*  86.4  * t>[NUM_CELLS  - 1]  * D0  * Vsample /SQR(L)); 
print/ ("\t7,6.  llf\n",  -loglO  (cation[H].c[lim  SURFACES  - 1])); 

} 

This  code  is  used  in  section  2. 


77.  Print  termination  information  to  include  the  day  of  termination  and  the  depth  of  chloride  and 
sulfate  penetration. 

(Print  termination  information  77)  = 
pnntf  ("\ n\n"); 

if  ( Day  > MaxDay  .value  A joint-change.flag  = FALSE)  prin//("ExceededuTIHEulimit . \n\n"); 
if  (Day  > MaxDay  .value  A joint.change.flag  = TRUE) 
pnn*/("ExceededuJOINTuLIFETIME . \n\nM); 
if  (sulfate.depth  > Max  Depth)  print/ ("Suit  ateuFailure . \n\n"); 
if  (chloride-depth  *L>  rebar.depth)  prinf/("ChlorideuFailure.\n\n"); 
pnn//("\"T\"\t'/,8 . 31f \n",  Time); 
pnntf  ("\" Day\"\t'/,8 . Olf \n",  Day); 
pnn^/("\"S04u(m)\"\t'/,8 . 31f  \n",  sulfate.depth  * L); 
pnn</("\"Clu(m)\"\t,/,8 . 311\n" , chlonde.depth  * L); 
pnntf  ("\ n\n"); 

This  code  is  used  in  section  2. 


78.  Print  final  system  state  for  any  desired  ion  or  salt. 

(Print  desired  output  78)  = 

pnntf  (''  FinaluSystemustate : \n\n" ); 

pnntf  ("\"L(m)\"u\t"); 

print/ ("\"  PsiV'uYt"); 

prjnf/("\"vD\"u\t"); 

print/  ("\"xi\"u\t"); 

print/ ( "\"phi\"\t"); 

pri7if/("\"pH\"u\t"); 

print/  ("V'f  c\"u\t"); 

for  (i  = 0;  i < num. cations;  i++) 

if  (cation[i\. output-flag  = TRUE)  pnntf  ("\"'/,s\" u\t'\  cation[i\.name); 
for  (j  = 0 ; j < num.anions;  j++) 

if  (anion  [j]. output-flag  = TRUE)  print/ ("\"7,s\"  u\t" , anion[j].name); 
for  (i  = 0;  i < num.cations;  i++) 
for  (j  = 0;  j < num.anions;  j++) 
if  (sol-array  [i][j]. output-flag  = TRUE) 

pnnf/("\"'/.s,/.s\"u\t",  cation[i].name , anion\j].name); 
pnntf  ("\  n"); 

for  (k  = 0;  k < NUM.SURFACES;  k++)  { 
print/ ("7,7 . 41i\t" , L * k); 
print/ ("7,7 . 31f  \t",  ’fffc]); 
if  (k  < NUM.CELLS)  { 

print/ ("7,6. 31iu\t",v[k]); 
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} 

else  { 

pnn//(,7.6.31fLj\t",  v[NUM_CELLS  - 1]); 

} 

print/  {"7,6 . 41f  u\t",  £'[£]); 

print/ {"7,6  AliuW' , </>'"[£]); 

print/ ("7,  6.31fu\t",  14.0  + logl  0 {anion[0].c[k]))\ 

print/  {''7,6 . 01fu\t",  strength[k]  * {log  {<j)'n[k])/ log  {<j>)))\ 

for  (z  = 0;  i < num.cations ; z++) 

if  {cation[i\. output-flag  = TRUE)  pmz//("'/,6.41fu\t",  ca/zozz  [z'].c[fc]); 
for  {j  = 0 ; j < num-anions ; j++) 

if  {anion  [j], output-flag  = TRUE)  print/ {''7,6 . 41f  u\t",  anion\j].c[k]); 
for  (z  = 0;  z < num.cations ; z++) 
for  (j  = 0 ; j < num.amons ; j++) 

if  (so/_arrap[z][/].ou/pzzL/?a<7  = TRUE)  pnnt/{"7,6 . 31f  u\t",  so/_arraj/[z][/].s[A:]); 
print/ {"\ n"); 

} 

elapsed-time  = biostime  {GET  _TIKE,  0 l)  — elapsed.time ; 
pnzi</(M\"Elapsedutimeu(sec)  : \"\t'/,6 . If \n",  elapsed.time /CLK_TCK); 

This  code  is  used  in  sections  2,  24,  and  26. 


